Add segment explaining best practices for CAD geometries
[dealii.wiki.git] / Frequently-Asked-Questions.md
1 # The deal.II FAQ
2 This page collects a few answers to questions that have frequently been asked about deal.II and that we thought are worth recording as they may be useful to others as well.
3
4 ## Table of Contents
5   * [The deal.II FAQ](#the-dealii-faq)
6     * [Table of Contents](#table-of-contents)
7     * [General questions on deal.II](#general-questions-on-dealii)
8       * [Can I use/implement triangles/tetrahedra in deal.II?](#can-i-useimplement-trianglestetrahedra-in-dealii)
9       * [I'm stuck!](#im-stuck)
10       * [I'm not sure the mailing list is the right place to ask ...](#im-not-sure-the-mailing-list-is-the-right-place-to-ask-)
11       * [How fast is deal.II?](#how-fast-is-dealii)
12       * [deal.II programs behave differently in 1d than in 2/3d](#dealii-programs-behave-differently-in-1d-than-in-23d)
13       * [I want to use deal.II for work in my company. Do I need a special license?](#i-want-to-use-dealii-for-work-in-my-company-do-i-need-a-special-license)
14     * [Supported System Architectures](#supported-system-architectures)
15       * [Can I use deal.II on a Windows platform?](#can-i-use-dealii-on-a-windows-platform)
16         * [Run deal.II in the Windows Subsystem for Linux](#run-dealii-in-the-windows-subsystem-for-linux)
17         * [Run deal.II natively on Windows](#run-dealii-natively-on-windows)
18         * [Run deal.II through a virtual box](#run-dealii-through-a-virtual-box)
19         * [Dual-boot your machine with Ubuntu](#dual-boot-your-machine-with-ubuntu)
20       * [Can I use deal.II on an Apple Macintosh?](#can-i-use-dealii-on-an-apple-macintosh)
21       * [Does deal.II support shared memory parallel computing?](#does-dealii-support-shared-memory-parallel-computing)
22       * [Does deal.II support parallel computing with message passing?](#does-dealii-support-parallel-computing-with-message-passing)
23       * [How does deal.II support multi-threading?](#how-does-dealii-support-multi-threading)
24       * [My deal.II installation links with the Threading Building Blocks (TBB) but doesn't appear to use multiple threads!](#my-dealii-installation-links-with-the-threading-building-blocks-tbb-but-doesnt-appear-to-use-multiple-threads)
25     * [Configuration and Compiling](#configuration-and-compiling)
26       * [Where do I start?](#where-do-i-start)
27       * [I tried to install deal.II on system X and it does not work](#i-tried-to-install-dealii-on-system-x-and-it-does-not-work)
28       * [How do I change the compiler?](#how-do-i-change-the-compiler)
29       * [I can configure and compile the library but installation fails. What is going on?](#i-can-configure-and-compile-the-library-but-installation-fails-what-is-going-on)
30       * [I get warnings during linking when compiling the library. What's wrong?](#i-get-warnings-during-linking-when-compiling-the-library-whats-wrong)
31       * [I can't seem to link/run with PETSc](#i-cant-seem-to-linkrun-with-petsc)
32         * [Is there a sure-fire way to compile deal.II with PETSc?](#is-there-a-sure-fire-way-to-compile-dealii-with-petsc)
33         * [I want to use HYPRE through PETSc](#i-want-to-use-hypre-through-petsc)
34         * [Is there a sure-fire way to compile dealii with SLEPc?](#is-there-a-sure-fire-way-to-compile-dealii-with-slepc)
35       * [Trilinos detection fails with an error in the file Sacado.hpp or <code>Sacado_cmath.hpp</code>](#trilinos-detection-fails-with-an-error-in-the-file-sacadohpp-or-sacado_cmathhpp)
36       * [My program links with some template parameters but not with others.](#my-program-links-with-some-template-parameters-but-not-with-others)
37       * [When trying to run my program on Mac OS X, I get image errors.](#when-trying-to-run-my-program-on-mac-os-x-i-get-image-errors)
38     * [C++ questions](#c-questions)
39       * [What integrated development environment (IDE) works well with deal.II?](#what-integrated-development-environment-ide-works-well-with-dealii)
40       * [Is there a good introduction to C++?](#is-there-a-good-introduction-to-c)
41       * [Are there features of C++ that you avoid in deal.II?](#are-there-features-of-c-that-you-avoid-in-dealii)
42       * [Why use templates for the space dimension?](#why-use-templates-for-the-space-dimension)
43       * [Doesn't it take forever to compile templates?](#doesnt-it-take-forever-to-compile-templates)
44       * [Why do I need to use typename in all these templates?](#why-do-i-need-to-use-typename-in-all-these-templates)
45       * [Why do I need to use this-&gt; in all these templates?](#why-do-i-need-to-use-this--in-all-these-templates)
46       * [Does deal.II require C++11 support?](#does-dealii-require-c11-support)
47         * [deal.II version 9.0.0](#dealii-version-900)
48         * [deal.II version 8.5.0 and previous](#dealii-version-850-and-previous)
49       * [Can I convert Triangulation cell iterators to DoFHandler cell iterators?](#can-i-convert-triangulation-cell-iterators-to-dofhandler-cell-iterators)
50     * [Questions about specific behavior of parts of deal.II](#questions-about-specific-behavior-of-parts-of-dealii)
51       * [How do I create the mesh for my problem?](#how-do-i-create-the-mesh-for-my-problem)
52       * [How do I describe complex boundaries?](#how-do-i-describe-complex-boundaries)
53         * [How do I ensure that my refined grid is correct when using CAD geometry?](#how-do-i-ensure-that-my-refined-grid-is-correct-when-using-cad-geometry)
54       * [How do I get the degree of freedom indices at vertices?](#how-do-i-get-the-degree-of-freedom-indices-at-vertices)
55       * [I am using discontinuous Lagrange elements (FE_DGQ) but they don't seem to have vertex degrees of freedom!?](#i-am-using-discontinuous-lagrange-elements-fe_dgq-but-they-dont-seem-to-have-vertex-degrees-of-freedom)
56       * [How do I access values of discontinuous elements at vertices?](#how-do-i-access-values-of-discontinuous-elements-at-vertices)
57       * [Does deal.II support anisotropic finite element shape functions?](#does-dealii-support-anisotropic-finite-element-shape-functions)
58       * [The graphical output files don't make sense to me -- they seem to have too many degrees of freedom!](#the-graphical-output-files-dont-make-sense-to-me----they-seem-to-have-too-many-degrees-of-freedom)
59       * [In my graphical output, the solution appears discontinuous at hanging nodes](#in-my-graphical-output-the-solution-appears-discontinuous-at-hanging-nodes)
60       * [When I run the tutorial programs, I get slightly different results](#when-i-run-the-tutorial-programs-i-get-slightly-different-results)
61       * [How do I access the whole vector in a parallel MPI computation?](#how-do-i-access-the-whole-vector-in-a-parallel-mpi-computation)
62       * [How to get the (mapped) position of support points of my element?](#how-to-get-the-mapped-position-of-support-points-of-my-element)
63     * [Debugging deal.II applications](#debugging-dealii-applications)
64       * [I don't have a whole lot of experience programming large-scale software. Any recommendations?](#i-dont-have-a-whole-lot-of-experience-programming-large-scale-software-any-recommendations)
65       * [Are there strategies to avoid bugs in the first place?](#are-there-strategies-to-avoid-bugs-in-the-first-place)
66       * [How can deal.II help me find bugs?](#how-can-dealii-help-me-find-bugs)
67       * [Should I use a debugger?](#should-i-use-a-debugger)
68       * [deal.II aborts my program with an error message](#dealii-aborts-my-program-with-an-error-message)
69       * [The program aborts saying that an exception was thrown, but I can't find out where](#the-program-aborts-saying-that-an-exception-was-thrown-but-i-cant-find-out-where)
70       * [I get an exception in virtual dealii::Subscriptor::~Subscriptor() that makes no sense to me!](#i-get-an-exception-in-virtual-dealiisubscriptorsubscriptor-that-makes-no-sense-to-me)
71       * [I get an error that the solver doesn't converge. But which solver?](#i-get-an-error-that-the-solver-doesnt-converge-but-which-solver)
72       * [How do I know whether my finite element solution is correct? (Or: What is the "Method of Manufactured Solutions"?)](#how-do-i-know-whether-my-finite-element-solution-is-correct-or-what-is-the-method-of-manufactured-solutions)
73       * [My program doesn't produce the expected output!](#my-program-doesnt-produce-the-expected-output)
74       * [The solution converges initially, but the error doesn't go down below 10<sup>-8</sup>!](#the-solution-converges-initially-but-the-error-doesnt-go-down-below-10-8)
75       * [My code converges with one version of deal.II but not with another](#my-code-converges-with-one-version-of-dealii-but-not-with-another)
76       * [My time dependent solver does not produce the correct answer!](#my-time-dependent-solver-does-not-produce-the-correct-answer)
77       * [My Newton method for a nonlinear problem does not converge (or converges too slowly)!](#my-newton-method-for-a-nonlinear-problem-does-not-converge-or-converges-too-slowly)
78       * [Printing deal.II data types in debuggers is barely readable!](#printing-dealii-data-types-in-debuggers-is-barely-readable)
79       * [My program is slow!](#my-program-is-slow)
80       * [How do I debug MPI programs?](#how-do-i-debug-mpi-programs)
81       * [I have an MPI program that hangs](#i-have-an-mpi-program-that-hangs)
82       * [One statement/block/function in my MPI program takes a long time](#one-statementblockfunction-in-my-mpi-program-takes-a-long-time)
83     * [I have a special kind of equation!](#i-have-a-special-kind-of-equation)
84       * [Where do I start?](#where-do-i-start-1)
85       * [Can I solve my particular problem?](#can-i-solve-my-particular-problem)
86       * [Why use deal.II instead of writing my application from scratch?](#why-use-dealii-instead-of-writing-my-application-from-scratch)
87       * [Can I solve problems over complex numbers?](#can-i-solve-problems-over-complex-numbers)
88       * [How can I solve a problem with a system of PDEs instead of a single equation?](#how-can-i-solve-a-problem-with-a-system-of-pdes-instead-of-a-single-equation)
89       * [Is it possible to use different models/equations on different parts of the domain?](#is-it-possible-to-use-different-modelsequations-on-different-parts-of-the-domain)
90       * [Can I solve problems with higher regularity requirements?](#can-i-solve-problems-with-higher-regularity-requirements)
91       * [Where do I start to implement a new Finite Element Class?](#where-do-i-start-to-implement-a-new-finite-element-class)
92     * [General finite element questions](#general-finite-element-questions)
93       * [How do I compute the error](#how-do-i-compute-the-error)
94       * [How to plot the error as a pointwise function](#how-to-plot-the-error-as-a-pointwise-function)
95       * [I'm trying to plot the right hand side vector but it doesn't seem to make sense!](#im-trying-to-plot-the-right-hand-side-vector-but-it-doesnt-seem-to-make-sense)
96       * [What does XXX mean?](#what-does-xxx-mean)
97     * [I want to contribute to the development of deal.II!](#i-want-to-contribute-to-the-development-of-dealii)
98     * [I found a typo or a bug and fixed it on my machine. How do I get it included in deal.II?](#i-found-a-typo-or-a-bug-and-fixed-it-on-my-machine-how-do-i-get-it-included-in-dealii)
99     * [I'm fluent in deal.II, are there jobs for me?](#im-fluent-in-dealii-are-there-jobs-for-me)
100
101 ## General questions on deal.II
102
103 ### Can I use/implement triangles/tetrahedra in deal.II?
104
105 This is truly one of the most frequently asked questions. The short answer
106 is: No, you can't. deal.II's basic data structures are too much tailored to
107 quadrilaterals and hexahedra to make this trivially possible. Implementing
108 other reference cells such as triangles and tetrahedra amounts to
109 re-implementing nearly all grid and DoF classes from scratch, along with
110 the finite element shape functions, mappings, quadratures and a whole host
111 of other things. Making triangles and tetrahedra work would certainly involve
112 having to write several ten thousand lines of code, and to make it usable
113 in all the rest of the library would require auditing a very significant
114 fraction of the 600,000 lines of code that make up deal.II today.
115
116 That said, the current specialization on quadrilaterals and hexahedra has
117 two very positive aspects: First, quadrilaterals and hexahedra typically
118 provide a significantly better approximation quality than triangular meshes
119 with the same number of degrees of freedom; you therefore get more accurate
120 solutions for the same amount of work. Secondly, because the shape of cells
121 are known, we can make a lot of things known to the compiler (such as the
122 number of iterations of a loop over all vertices of a cell) which avoid a
123 large number of run-time computations and makes the library as fast as it
124 is. A simple example is that in deal.II we know that a loop over all
125 vertices of a cell has exactly `GeometryInfo<dim>::vertices_per_cell`
126 iterations, a number that is known to the compiler at compile-time. If we
127 allowed both triangles and quadrilaterals, the loop would have
128 `cell->n_vertices()` iterations, but this would in general not be known at
129 compile time and consequently not allow the compiler to optimize on.
130
131 If you do need to work with a geometry for which all you have is a
132 triangular or tetrahedral mesh, then you can convert this mesh into one
133 that consists of quadrilaterals and hexahedra using the `tethex` program,
134 see https://github.com/martemyev/tethex .
135
136 ### I'm stuck!
137
138 Further down below on this page (in the debugging section) we list a number
139 of strategies on how to find errors in your program. If your question is
140 how to implement something new for which you don't know where to start,
141 have you taken a look at the set of tutorial programs and checked whether
142 one or the other already has something that's close to what you want?
143
144 That said, there will be situations where documentation doesn't help and
145 where you need other someone else's opinion. That's what the [deal.II
146 mailing lists](http://dealii.org/mail.html) are there for: Feel free to
147 ask! You may also wish to subscribe to the users' list -- not so much
148 because someone else might ask the same question you have, but because
149 reading the list gives you background information on things others are
150 working on that may help you when you want to do something similar.
151
152 When asking for help on the mailing list, be specific. We frequently get mail of the following kind:
153 <pre>
154   I'm trying to do X. This works fine but it fails when I try to transfer
155   the data to my MyClass::Estimator object. I tried to use something
156   similar to what's done in a couple of tutorial programs but it doesn't
157   work. I'm new at C++ and I just can't seem to get the syntax right.
158 </pre>
159
160 This message doesn't contain nearly enough information for anyone to really
161 help you: we don't know what `MyClass::Estimator` is, we don't know how you
162 try to transfer data, we haven't seen your code, and we haven't seen the
163 compiler's error messages. (For more examples of how not to write help
164 requests, see [Section 3.2 of this
165 document](http://faculty.washington.edu/dchinn/how-not-to-code.pdf).) We
166 could poke in the dark, but it would probably be more productive if you
167 gave us a bit more detail explaining what doesn't work: show us the code
168 you implemented, show us the compiler's error message, or be specific in
169 some other way in describing what the problem is!
170
171 ### I'm not sure the mailing list is the right place to ask ...
172
173 Yes, it probably is. Please direct your questions to the mailing list and
174 not to individual developers. There are many reasons:
175
176   1. Others might have similar questions in the future and can search the
177      archives.
178   1. There are many active users on the mailing list that are happy to
179      help. There probably is someone who did something very similar before.
180   1. Imagine everyone would stop using the mailing list and email us
181      directly. We would spend most of our time answering the same questions
182      over and over.
183   1. Many users are reading the mailing list and are interested in deal.II
184      in general and are learning by skimming emails. Give them a chance.
185   1. As a consequence of all this, we typically prioritize questions on
186      mailing lists over emails sent directly to us asking for help.
187   1. Don't be afraid. There are no stupid questions (only off-topic ones).
188      Everyone started out at some point. Asking the questions in the open
189      helps us improve the library and documentation.
190
191 That said, if there is something you can not discuss in the open, feel free
192 to contact us!
193
194
195 ### How fast is deal.II?
196
197 The answer to this question really depends on your metric. If you had to
198 write, say, a Stokes solver with a particular linear solver, a particular
199 time stepping scheme, on a piecewise polygonal domain, and Q2/Q1 elements,
200 you can write a code that is 20% or 30% faster than what you would get when
201 using deal.II because you know the building blocks, shape functions,
202 mappings, etc. But it'll take you 6 months to do so, and 20,000 lines of
203 code. On the other hand, when using deal.II, you can do it in 2 weeks and
204 204 lines (that's the number of semicolons in step-22).
205
206 In other words, if by "fast" you mean the absolute maximal efficiency in
207 terms of CPU time deal.II is more than likely to lose against a
208 hand-written Fortran77 code. But for most of us, the real question of
209 "fast" also includes the time it takes to get the code running and
210 verified, and in that case deal.II is most likely the fastest library out
211 there simply by virtue of the fact that it is by far the largest and most
212 comprehensive finite element library available as Open Source.
213
214 This all, by the way, does not mean that we don't care about speed: We
215 spend a lot of effort profiling the library and working on the hot spots to
216 make codes fast. The discussion of this issue in the introduction of
217 step-22 is a good example. There are also some guidelines below on how to
218 profile your code in the debugging section of this FAQ.
219
220 ### deal.II programs behave differently in 1d than in 2/3d
221
222 In deal.II, you can write programs that look exactly the same in 2d and 3d,
223 but there are cases where 1d is slightly different. That said, this is an
224 area that we have significantly rewritten, and starting with deal.II 7.1,
225 most cases should work in 1d in just the same way as they do in 2d/3d. If
226 you find something that doesn't work, please report it to the mailing list.
227
228 Historically, the differences primarily resulted from the fact that in
229 deal.II, we represent vertices differently from lines and quads; whereas
230 the latter can store information (for example boundary indicators, user
231 flags, etc) vertices don't. As a consequence, the boundary indicator of a
232 boundary part in 1d (i.e. either the left or right vertex) were determined
233 by convention, rather than by setting it explicitly: the left boundary of a
234 1d domain always had boundary indicator zero, the right boundary always
235 boundary indicator one. This was different from the 2d/3d case where by
236 default (unless you explicitly set things differently) all boundaries have
237 indicator zero. This left-boundary-has-id-0, right-boundary-has-id-1 is
238 still the default today, but at least you can set the boundary indicators
239 of these end-points to something different today.
240
241 A second difference is that vertices have no extent, and so you can't apply
242 quadrature to them. As a consequence, the FEFaceValues class wasn't usable
243 in 1d. Again, this should work these days: every quadrature formula that
244 has a single quadrature point is a valid one for points as well.
245
246 ### I want to use deal.II for work in my company. Do I need a special license?
247
248 Before going into any more details, you **need** to carefully read the
249 license deal.II is under. In particular, the explanations below are not meant
250 to be legal advice and does not override the provisions in the Open Source
251 license.
252
253 However, before this, let us provide our overarching philosophy: It is our
254 intention to have constructive relationships with those who want to use our
255 work commercially, and we encourage commercial use. After having used a
256 more restrictive license until 2013, we have come to the conclusion that
257 these licenses serve neither side particularly well: it made commercial use
258 difficult, and the lack of commercial use deprived us of critical feedback,
259 potential contributions from professional users, and our users of potential
260 employment opportunities. Everyone is better off with the LGPL license we
261 are using now, and we hope that deal.II also finds use in commercial
262 settings.
263
264 Now for the smaller print: Generally, the LGPL is a fairly liberal license.
265 In particular, if you *develop a code based on deal.II*, then there is no
266 requirement that you also open source your own code: you can keep it closed
267 source, under a proprietary license, and you don't need to give it to
268 anyone (neither your customers nor to us).
269
270 The LGPL is only restrictive in that the *changes you make to deal.II
271 itself* must also be licensed under the LGPL. There is not frequently a
272 need to change the library itself, and in many of these cases you will
273 probably be interested to get them into the upstream development sources
274 anyway (e.g., in cases of bugs) rather than having to forward port them
275 indefinitely. Of course, we are interested in this as well. However, there
276 is no such requirement that you upstream these changes: the only people you
277 have to make these modifications to deal.II available to are your
278 customers.
279
280 As mentioned above, the preceding paragraphs are not a legal
281 interpretation. For definite interpretations of the LGPL, you may want to
282 consult lawyers familiar with the topic or search the web for more detailed
283 interpretations.
284
285
286 ## Supported System Architectures
287
288 ### Can I use deal.II on a Windows platform?
289
290 deal.II has been developed with a Unix-like environment in mind and it
291 shows in a number of places regarding the build system and compilers
292 supported. That said, there are multiple methods to get deal.II running if
293 you have a Windows machine.
294
295 #### Run deal.II in the Windows Subsystem for Linux
296
297 Windows 10 has gained a compatibility layer for running Linux binaries
298 natively on Windows. You can find more information on the
299 [Wikipedia page](https://en.wikipedia.org/wiki/Windows_Subsystem_for_Linux).
300 This means you do not have to use
301 [virtualization](#run-dealii-through-a-virtual-box), or [dual
302 boot](#dual-boot-your-machine-with-ubuntu) any more to install a
303 full-featured Linux distribution! We summarize the installation on a
304 separate wiki page on [[Windows]].
305
306 #### Run deal.II natively on Windows
307
308 Since deal.II 8.4.0 we have experimental support for Microsoft Visual Studio (2013 and 2015). See the separate page on [[Windows]] for more details.
309
310 #### Run deal.II through a virtual box
311
312 The simplest way to try out deal.II is to run it in a premade virtual
313 machine. You can download the virtual machine for VirtualBox from
314 http://www.math.clemson.edu/~heister/dealvm/ and run it inside windows.
315
316 Note that your experience depends on how powerful your machine is. More
317 than 4GB RAM are recommended. A native installation of Linux is preferable
318 (see below).
319
320 #### Dual-boot your machine with Ubuntu
321
322 The simplest way to install Linux as a Windows user is to dual-boot.
323 Dual-boot means that you simply install a second operating system on your
324 computer and you choose which one to start when you boot the machine. Most
325 versions of Linux support installing themselves as a second operating
326 system. One example is using the Ubuntu installer for Windows. This
327 installer will automatically dual-boot your system for you in a safe and
328 fully reversible manner. Simply follow the instructions on
329 http://www.ubuntu.com/download/desktop/install-ubuntu-with-windows
330
331 If at some point in the future you wish to remove Ubuntu from your system,
332 from the Windows program manager (add-remove programs in older versions and
333 programs and features in newer versions) you can simply uninstall Ubuntu as
334 you would any other program.
335
336 *Note:* The actual install file is linked through the text "Windows
337 installer" in the first gray box.  You will be prompted to donate to
338 Ubuntu, which is entirely optional. You will also be prompted to use a
339 different version of Ubuntu if you use Windows 8.
340
341 ### Can I use deal.II on an Apple Macintosh?
342
343 Yes, at least on the more modern OS X operating systems this works just
344 fine. deal.II supports native compilers shipping with XCode as well as gcc
345 from Mac Ports.
346
347 The only issue we are currently aware of is that if deal.II is configured
348 to interface with PETSc, then PETSc needs to be configured with the
349 <code>--with-x=0</code> flag to prevent linking in the X11 libraries (you
350 probably won't need them anyway). Installing with PETSc has a myriad of
351 other problems, though we believe that we have a way to stably interface
352 it. You may want to read through the PETSc-related entries further down,
353 however.
354
355 ### Does deal.II support shared memory parallel computing?
356
357 Yes. deal.II supports multithreading with the help of the
358 [http://www.threadingbuildingblocks.org Threading Building Blocks (TBB)
359 library](c967ec2ff74d85bd4327f9f773a93af3]). It is enabled by default and
360 can be controlled via the `DEAL_II_WITH_THREADS` configuration toggle
361 passed to `cmake` (see the deal.II readme file).
362
363 ### Does deal.II support parallel computing with message passing?
364
365 Yes, and in fact it has been shown to scale very nicely to at least 16,384
366 processor cores in a paper by Bangerth, Burstedde, Heister and Kronbichler.
367 You should take a look at the documentation modules discussing parallel
368 computing, as well as the step-40 tutorial program.
369
370
371 ### How does deal.II support multi-threading?
372
373 deal.II will use multi-threading using several approaches:
374 1. some BLAS routines might be multi-threaded (typically using OpenMP).
375    This can be controlled from the command line using OMP_NUM_THREADS (also
376    see the entry in the FAQ below)
377 2. Many places in the library are parallelized using the Threading Building
378    Blocks (TBB) library.
379
380 MPI_InitFinalize() has an optional third argument that specifies the number
381 of threads to use for the TBB. The default is 1. This gets send to the TBB
382 via a call to  MultithreadInfo::set_thread_limit(). If you pass
383 numbers::invalid_unsigned_int into MPI_InitFinalize (or if you don't use
384 that class, call set_thread_limit directly) then TBB will use the maximum
385 number of threads that makes sense (and you can limit it using
386 DEAL_II_NUM_THREADS from the command line).
387
388 Also note that while our Trilinos wrappers support multi-threading, the
389 PETSc wrappers do not support this at this time, so you need to run with
390 one thread per process.
391
392 ### My deal.II installation links with the Threading Building Blocks (TBB) but doesn't appear to use multiple threads!
393
394 This may be a quirky interaction with the [GOTO
395 BLAS](http://www.tacc.utexas.edu/tacc-projects/gotoblas2/) :-( If you use
396 Trilinos or PETSc, both of these require a BLAS library from your system,
397 and the deal.II cmake configuration will make sure that it is linked with.
398 The problem stems from the fact that by default, the GOTO BLAS will simply
399 grab all cores of the system for its own use, and -- before your `main()`
400 function even starts, allow the main thread to use only a single core. (For
401 the technically interested: it sets the processor scheduling affinity mask,
402 using `set_sched_affinity` to a single bit.)
403
404 When the TBB initialization runs, still before `main()` starts, it will
405 find that it can only run on a single core and will consequently not be
406 able to work on multiple tasks in parallel.
407
408 The solution to this problem is to forbid the GOTO BLAS to grab all
409 processors for itself, since we spend very little time in BLAS anyway. This
410 can be done by setting either the `OMP_NUM_THREADS` or `GOTO_NUM_THREADS`
411 environment variables to 1, see
412 http://www.tacc.utexas.edu/tacc-software/gotoblas2/faq .
413
414
415 ## Configuration and Compiling
416
417 ### Where do I start?
418
419 Have a look at the  [ReadMe instructions](http://www.dealii.org/developer/readme.html) for details on how to configure and install the library with `cmake`.
420
421 ### I tried to install deal.II on system X and it does not work
422
423 That does occasionally (though relatively rarely) happen, in particular if
424 you work on an operating system or with a compiler that the primary
425 developers don't have access to. In a case like this, you should ask for
426 help on the mailing list. However, remember: If your question only contains
427 the text "I tried to install deal.II on system X and it does not work" then
428 that's not quite enough to figure out what is happening. Even though the
429 people developing this software belong to the most able programmers in the
430 universe (and a decent number of parallel universes), all of us need data
431 to find errors. So, whatever went wrong, paste the error message into your
432 email. If the error is from the `cmake` invocation, show us the error
433 message that was printed on screen.
434 If the error happens after configuring and during compiling, add lines from
435 screen output showing the error to the mail.
436
437
438 ### How do I change the compiler?
439
440 deal.II can be compiled by a number of compilers without problems (see the
441 section [prerequisites](http://www.dealii.org/readme.html#prerequisites) in
442 the readme file). If `cmake` does not pick the right one, selecting another
443 is simple, and described in a
444 [section](http://www.dealii.org/developer/development/cmake.html#compiler)
445 in the [cmake
446 documentation](http://www.dealii.org/developer/development/cmake.html).
447
448 ### I can configure and compile the library but installation fails. What is going on?
449
450 If you configure with the default ``CMAKE_INSTALL_PREFIX``, the library is configured to installed to ``/usr/local`` and this fails without superuser rights with an error message like
451 ```
452 CMake Error at cmake/scripts/cmake_install.cmake:42 (FILE):
453   file cannot create directory: /usr/local/common/scripts.  Maybe need
454   administrative privileges.
455 ```
456 Please see the [readme](http://www.dealii.org/developer/readme.html#configuration) on how to pick an install directory with write access (for example some path below your home directory).
457
458 ### I get warnings during linking when compiling the library. What's wrong?
459
460 On some linux distributions with particular versions of the system
461 compiler, one can get warnings like these during the linking stage of
462 compiling the library:
463 ```
464 `.L3019' referenced in section `.rodata' of
465 /home/bangerth/deal.II/lib/lac/sparse_matrix.float.g.o: defined in discarded section
466 `.gnu.linkonce.t._ZN15SparsityPattern21optimized_lower_boundEPKjS1_RS0_'
467 of /home/bangerth/deal.II/lib/lac/sparse_matrix.float.g.o
468 ```
469
470 While annoying, these warnings do not actually seem to indicate anything
471 particularly harmful. Apparently, the compiler generates the same code
472 multiple times in exactly the same form, and the linker is only warning
473 that it is throwing away all but one of the copies. There doesn't seem to
474 be way to avoid these warnings, but they can be safely ignored.
475
476 ### I can't seem to link/run with PETSc
477
478 Recent deal.II releases support PETSc 3.0 and later. This works, but there
479 are a number of things that can go wrong and that result in compilation or
480 linker errors, as explained below. If your program links properly with
481 PETSc support, it will very likely also produce the correct results.
482
483 If you get errors like this when trying to run step-17 of the tutorials,
484 even though linking seems to have succeeded just fine:
485 ```
486    [make run
487    ============================ Running step-17
488    ./step-17: error while loading shared libraries: libpetsc.so: cannot open
489               shared object file: No such file or directory
490    make: *** [run](step-17]) Error 127
491 ```
492
493 this means is that while linking, the compiler could find the libpetsc.so
494 library, but the executable can't find it when running. The reason is that
495 we can tell the linker where to look, but the executable apparently did not
496 remember this (this is the standard Unix behavior). What you have to do is
497 to set the LD_LIBRARY_PATH to include the path to the PETSc libraries. For
498 example, under `bash` you would have to do this:
499 ```
500    export LD_LIBRARY_PATH=/path/to/petsc/libraries:$LD_LIBRARY_PATH
501 ```
502
503 If you do so, the Unix loader can query the environment variable for where
504 to find this particular library when trying to run the executable, and
505 running the program should succeed.
506
507 Similarly, if you get errors of the kind during linking
508 ```
509 /home/xxx/deal.II/lib/libdeal_II.g.so: undefined reference to
510 `KSPSetInitialGuessNonzero(_p_KSP*, PetscTruth)'
511 /home/xxx/deal.II/lib/libdeal_II.g.so: undefined reference to
512 `VecAXPY(_p_Vec*, double, _p_Vec*)'
513 ...
514 ```
515
516 then the compiler can't seem to find the PETSc libraries. The solution is
517 as above: specify the path to those libraries via `LD_LIBRARY_PATH`.
518
519
520 #### Is there a sure-fire way to compile deal.II with PETSc?
521
522 Short answer is "No". The slightly longer answer is, "PETSc has too many
523 knobs, switches, dials, and a kitchen sink too many for its own damned
524 good. There is not a sure-fire way to compile deal.II with PETSc!". It
525 turns out that PETSc is a very versatile machine and, as such, there is no
526 shortage of things that can go wrong in trying to configure PETSc to work
527 seamlessly with deal.II on a first attempt. We have all struggled with
528 this, although it has become a lot better in recent years.
529
530 You can find instructions on how to install PETSc linked to from the
531 deal.II ReadMe file, or going directly to
532 http://www.dealii.org/developer/external-libs/petsc.html .
533
534 #### I want to use HYPRE through PETSc
535
536 Hypre implements algebraic multigrid methods (AMG) as preconditioners, for
537 example the BoomerAMG method. AMGs are among the most efficient
538 preconditioners available and they have also been shown to be scalable to
539 thousands of processors. deal.II allows the use of Hypre through the
540 PETScWrappers::PreconditionBoomerAMG class; it is used in `step-40`. Hypre
541 can be installed as a sub-package of PETSc and deal.II can access it
542 through the PETSc interfaces.
543
544 To use the Hypre interfaces through PETSc, you need to configure PETSc as
545 discussed in http://www.dealii.org/developer/external-libs/petsc.html  ,
546 and add the following switch to the command line: `--download-hypre=1`.
547
548 #### Is there a sure-fire way to compile dealii with SLEPc?
549
550 Happily, the answer to this question is a definite yes; that is, <b>if you
551 have successfully compiled and linked PETSc already</b>.
552
553 The real trick here is that during configuration SLEPc will pull out
554 PETSc's configuration and just does whatever that tells it to do. Detailed
555 steps are discussed in
556 http://www.dealii.org/developer/external-libs/slepc.html .
557
558 Once deal.II is compiled, it is worth to start by looking at the step-36
559 tutorial program to see how to get started using the interface with SLEPc.
560
561 <i>
562 Note: To use the solvers and other algorithms SLEPc provides it is
563 absolutely essential to have your PETSc installation working correctly
564 since they share the same vector-matrix (and other) data structures.
565 </i>
566
567
568 ### Trilinos detection fails with an error in the file `Sacado.hpp` or `Sacado_cmath.hpp`
569
570 This is a complicated one (and it should also be fixed in more recent
571 Trilinos versions). In the Trilinos file `Sacado_cmath.hpp`, there is some
572 code of the form
573 ```cpp
574   namespace std
575   {
576     inline float acosh(float x)
577     {
578       return std::log(x + std::sqrt(x*x - float(1.0)));
579     }
580     ...
581   }
582 ```
583
584 In other words, Sacado is putting things into namespace `std`. The functions it
585 is putting there are functions that have been defined by the C99 standard but
586 that didn't make it into the C++98 standard before; some of them are widely
587 used. The problem is that these functions were later added to the standard and
588 so if your compiler is new enough (e.g. GCC 4.5 and later) then the compiler's
589 C++ standard library already contains these functions. Adding them again in this
590 file then yields errors of the kind
591 ```
592 /home/.../trilinos-10.4.2/include/Sacado_cmath.hpp: In function 'float std::acosh(float)':
593 /home/.../trilinos-10.4.2/include/Sacado_cmath.hpp:41:16: error: redefinition of 'float std::acosh(float)'
594 /usr/include/c++/4.5/tr1_impl/cmath:321:3: error: 'float std::acosh(float)' previously defined here
595 ```
596
597 The only useful way to avoid this error is to edit the Trilinos header
598 file. To do this, find and open the file `include/Sacado_cmath.hpp` in the
599 directory in which Trilinos was installed. Then change the block enclosed
600 in
601 ```cpp
602   namespace std
603   {
604     ...
605   }
606 ```
607
608 to read
609 ```cpp
610 #ifndef _GLIBCXX_USE_C99_MATH
611   namespace std
612   {
613     ...
614   }
615 #endif
616 ```
617
618 What this will do is make sure that the new members of namespace `std` are
619 only added if the compiler has not already done so itself.
620
621
622
623
624 ### My program links with some template parameters but not with others.
625
626 deal.II has many types for whose initialization you need to provide a
627 template parameter, e.g. `SparseMatrix<double>`. The implementation of
628 these classes can typically be found in files ending `.templates.h`, e.g.
629 `sparse_matrix.templates.h`. The corresponding `.cc` files, e.g.
630 `sparse_matrix.cc`, essentially only provide the explicit instantiations of
631 these classes for the most commonly used template parameters. Sometimes
632 this is done by including a corresponding `.inst` file, e.g.
633 `sparse_matrix.inst`.
634
635 If you want to use a data type with a template parameter for which there is
636 an explicit instantiation, you only need to include the respective `.h`
637 header file, e.g. `sparse_matrix.h`. If, however, you want to use a
638 template parameter for which there is no explicit instantiation in the
639 corresponding `.cc` file, you have to include the respective `.templates.h`
640 file in order for your program to link successfully.
641
642 The reason for all of this is essentially a matter of reducing compilation
643 time. As long as you use data types with template parameters for which
644 there is an explicit instantiation - and this should be the case most of
645 the time - you do not need to compile the respective (lengthy) .templates.h
646 file every time you compile your code. If, however, you need to use an
647 instance of e.g. `SparseMatrix<bool>`, you have to include the respective
648 `.templates.h` file and you have to compile it along with the remaining
649 files of your program every time.
650
651 ### When trying to run my program on Mac OS X, I get image errors.
652
653 You may encounter an error of the form
654
655 ```
656 dyld: Library not loaded: libdeal_II.g.7.0.0.dylib. Reason: image not found
657 ```
658
659 on OS X. This goes hand in hand with the following message you should have
660 gotten at the end of the output of `./configure`:
661
662 ```
663      Please add the line
664         export DYLD_LIBRARY_PATH=\$DYLD_LIBRARY_PATH:$DEAL2_DIR/lib
665      to your .bash_profile file so that OSX will be
666      able to find the deal.II shared libraries when
667      executing your programs.
668 ```
669
670 What happens is this: when you say "make all", all the deal.II files are
671 compiled and linked into a library (called libdeal_II.g.7.0.0) which on Macs
672 have the file ending .dylib. Then you go to examples/step-1 and compile your
673 program, which uses all the functions and classes that have previously been
674 put into this library.
675
676 Now the following happens: On most operating systems, the actual executable
677 program (i.e. the file step-1 in your directory that resulted from compiling)
678 does not contain any information that would indicate where the various
679 libraries that it uses can be found. For example, the step-1 program does not
680 know where the libdeal_II.g.7.0.0.dylib file is. This is just how most
681 operating systems function. But when you want to execute the program, somehow
682 the program has to know where the library it needs is located. On most
683 unix-like operating systems, this is done by setting an "environment
684 variable" -- on linux this would the variable "LD_LIBRARY_PATH", on Mac OS X
685 it is "DYLD_LIBRARY_PATH".
686
687 So to let the operating system know where the library is located, you could
688 type
689   export DYLD_LIBRARY_PATH=$DYLD_LIBRARY_PATH:/Users/renjun/deal.ii/lib
690 every time before you want to execute the program. That would be cumbersome. A
691 simpler way would be if this export command is executed every time when you
692 open a new shell window. This can be achieved by putting this command in a
693 file that is executed every time you open a shell window. Depending on what
694 shell you use, these files are alternatively called
695   .cshrc
696   .bashrc
697   .bash_profile
698 or similar. I'm not quite sure which file is relevant for you, but you can try
699 them one after the other by putting the text in there, closing the window,
700 opening it again, and then trying to execute
701   ./step-1
702 (or saying "make run" in this directory) and seeing whether that works.
703
704 ## C++ questions
705
706 ### What integrated development environment (IDE) works well with deal.II?
707
708 The short answer is probably: whatever works best for you. deal.II uses the
709 build tool CMake, which can generate a project description for virtually every
710 IDE. In the past, many of the main developers have used emacs (or even vi), but
711 there are much better tools around today, such as [eclipse](http://www.eclipse.org/),
712 [KDevelop](http://www.kdevelop.org),
713 [Xcode](http://developer.apple.com/technologies/tools/),
714 [QtCreator](http://qt.nokia.com/products/developer-tools/), all of which
715 have been used by people using deal.II.
716
717 We have gathered some notes on using the following IDEs for deal.II:
718   - [[Eclipse]]
719   - [[KDevelop]]
720   - [[emacs]]: While we don't recommend using emacs any more, this link provides a couple of notes on formatting styles used within deal.II.
721
722 When thinking about what IDE to use, keep this in mind: Many of us have
723 used emacs (or, worse, vi) for years and feel very comfortable with it.
724 But, emacs and vi were both started in 1976, at a time when computers had
725 little memory, virtually no CPU power, and only text-based interfaces.
726 While they have of course become a lot better over time, the design
727 limitations this involved are still very much part of the code base:
728 fundamentally, they are both still text-based and file-oriented. What IDEs
729 can provide are multiple views of the same project in graphical and textual
730 form and, more importantly, can integrate entire projects spanning hundreds
731 of files in multiple directories: they know where a variable is declared
732 (even if it's in a different file), what it's type is, and the properties
733 of this type. Neither emacs nor vi nor any other older editor can provide
734 anything that comes even close to what kdevelop or eclipse can offer in
735 this regard.
736
737 What all this implies is that you should consider using one of the more
738 modern tools, even if you're well acquainted with an existing, older one.
739 Of course it takes a while to get used to a new application but my
740 (Wolfgang's) experience with switching from emacs to kdevelop was that I
741 have become '''so''' much more productive by using modern tools that the
742 time invested in learning it was amortized very quickly. I found this
743 experience a real eye-opener!
744
745 ### Is there a good introduction to C++?
746
747 There are of course many good books and online resources that explain C++.
748 As far as websites are concerned,
749 [www.cplusplus.com](http://www.cplusplus.com) has both [reference material
750 for individual classes of the C++ standard
751 library](http://cplusplus.com/reference/) as well as a [a tutorial on parts
752 of the C++ language](http://cplusplus.com/tutorial) if you want to brush up
753 on the correct syntax of things.
754
755 ### Are there features of C++ that you avoid in deal.II?
756
757 There are few things that we avoid <i>as a matter of principle.</i> C++ is,
758 by and large, a pretty well designed language in the sense that its
759 features are there because they have been found to be useful by a lot of
760 people. As an example, people have found that it is easier to write and
761 debug code that throws exceptions in error cases rather than encoding error
762 situations by special return values (e.g. by returning -1). There are of
763 course ways to avoid exceptions (or templates, or certain parts of the C++
764 standard libraries, or any number of other things people have found
765 objectionable in C++) and some software projects have chosen to restrict
766 the use of C++ (for example Mozilla) or to emulate only those parts of C++
767 they like in C (e.g. the GNOME desktop environment, which leads to awkward
768 to understand code
769 [as described here](http://developer.gnome.org/gobject/stable/howto-gobject-methods.html)).
770
771 But ultimately, it is our belief that these approaches shoot their
772 inventors in the foot: they avoid features of C++ that were really intended
773 to make programming life simpler. It may be simpler for novice programmers
774 to read code without templates; ultimately, however, learning to read and
775 use templates will make you a much more productive programmer since you
776 don't write the same code multiple times. As a consequence, the use of C++
777 is driven by the question of what is best suited to write a particular
778 algorithm, not by abstract considerations. This fits into the realization
779 that deal.II is a large piece of software -- not a small research project
780 -- that requires professional software management practices and for which
781 long term development can no longer be driven by an individual programmer's
782 preferences of style.
783
784 ### Why use templates for the space dimension?
785
786 The fundamental motivation for this is to use dimension-independent
787 programming, i.e. you want to write code in such a way that it looks
788 exactly the same in 2d as in 3d (or 1d, for that matter). There are of
789 course many ways to do this (and libraries have done this for a long time
790 before deal.II has). The three most popular ones are to use a preprocessor
791 `#define` that sets the space dimension globally, to use a global variable
792 that does this, and to have each object have a member variable that denotes
793 the space dimension it is supposed to live in (in much the same way as the
794 template argument does in deal.II). Neither of these approaches is optimal
795 (nor is our own approach to use templates), however. In particular, using a
796 preprocessor symbol or a global variable will not allow you to mix and
797 match objects of different dimensionality. There are situations when you
798 want to do that; for example deal.II internally builds higher dimensional
799 quadrature formulas as tensor products of lower dimensional ones, and in
800 application codes you may wish to discretize both volume models (e.g.
801 simulating 3d models of plate tectonics and mountain belt formation) with
802 surface models (e.g. erosion processes on the 2d earth surface).
803
804 This leaves the option to have a member variable denoting the space
805 dimension in each object, a choice most other finite element libraries have
806 followed. But this isn't optimal either, for two reasons. For example,
807 consider this code that describes the equivalent of the `Point<dim>` class
808 for points in dim-dimensional space and its `norm()` member function:
809
810 ```cpp
811 class Point
812 {
813 public:
814   Point (const unsigned int dimension)
815     : dim(dimension),
816       coordinates (new double[dim])
817   {}
818
819   ~Point() { delete[] coordinates; }
820
821   double norm () const;
822   ...
823 private:
824   unsigned int dim;
825   double *coordinates;
826 };
827
828 double Point::norm () const
829 {
830   double s = 0;
831   for (unsigned int d=0; d<dim; ++d)
832     s += coordinates[d] * coordinates[d];
833   return std::sqrt(s);
834 }
835 ```
836
837 This is going to lead to rather slow code, for multiple reasons:
838
839  - The constructor and destructor have to allocate and deallocate memory on
840    the heap, both expensive processes.
841
842  - When accessing any element of the `coordinates` array, two pointers have
843    to be dereferenced. For example, the access to `coordinates[d]` really
844    expands to `*(this->coordinates + d)`.
845
846  - The compiler can not optimize the loop since the upper bound `dim` of
847    the loop variable is unknown at compile time.
848
849
850 Compare this to the way deal.II (approximately) implements this class:
851
852 ```cpp
853   template <int dim>
854   class Point
855   {
856     public:
857       Point () {}
858       ~Point() {}
859
860       double norm () const;
861       ...
862     private:
863       double coordinates[dim];
864   };
865
866   template <int dim>
867   double Point<dim>::norm () const
868   {
869     double s = 0;
870     for (unsigned int d=0; d<dim; ++d)
871       s += coordinates[d] * coordinates[d];
872     return std::sqrt(s);
873   }
874 ```
875
876 Here, the following holds:
877
878  - Constructor and destructor do not have to allocate and deallocate memory
879    on the heap; rather, since the size of the `coordinates` array is known
880    at compile time (i.e. whenever you instantiate the template for a
881    particular dimension), the array lives on the stack. It is also much
882    smaller than before: the dimension is encoded in the type and doesn't
883    need a memory location, we don't need to store a pointer to an array,
884    and we don't incur the memory overhead of having to manage an object on
885    the heap.
886
887  - When accessing any element of the `coordinates` array, only one pointer
888    has to be dereferenced. For example, the access to `coordinates` really
889    expands to `*(this + d)`.
890
891  - The compiler can optimize the loop since the upper bound `dim` of the
892    loop variable is known at compile time. In particular, for a point in
893    2d, the code the compiler will produce is likely to look more like this
894    because the loop can be unrolled and the loop counter can be optimized
895    away:
896 ```cpp
897 double Point<2>::norm () const
898 {
899   return std::sqrt(coordinates[0] * coordinates[0] + coordinates[1] * coordinates[1]);
900 }
901 ```
902 Obviously, for a 3d point, the code will look differently, but the compiler
903 can do this since it knows what the dimension of the point is at compile
904 time.
905
906 There is another reason for the deal.II way: type safety. In short, a 2d
907 point is not the same as a 3d point. If you assign one to the other, then
908 this may be on purpose and the executable should simply change the value of
909 the `dim` member variable from 2 to 3. But it may also be a legitimate
910 error -- for example, you shouldn't be able to use 2d points to initialize
911 the 3d quadrature points needed to integrate on a 3d cell. This can of
912 course be caught by run-time checks, but the reason for strongly typed
913 languages such as C++ has always been that it is much more efficient if the
914 compiler can already catch this sort of error at compile time. Using
915 templates for the space dimension avoids these sort of mistakes up front by
916 forcing the programmer to explicitly specify her intent, rather than
917 encoding intent in assertions.
918
919 Of course there are also downsides to using templates. Most notably, error
920 messages that involve templates are notoriously unreadable, and that
921 compiling template heavy code is slow: for example, we have to compile the
922 `Point` class three times (for dim=1, dim=2 and dim=3) rather than only
923 once. Nevertheless, we believe that these valid objections do not outweigh
924 the benefits of templates.
925
926 ### Doesn't it take forever to compile templates?
927
928 Yes, in general it does. The reason is that while for non-templates it is
929 enough to put the ''declaration'' of a function into the header file and
930 the ''definition'' into the `.cc` file, for templates that doesn't work.
931 Let's say you have something like
932 ```cpp
933   template <typename T> T square (const T & t);
934 ```
935
936 in your header file and you put the definition
937 ```cpp
938   template <typename T> T square (const T & t) { return t*t; }
939 ```
940
941 into the `.cc` file, then the compiler will say "Yes, I saw this template,
942 and if I see a use of this function later on I will generate a function
943 from it by replacing `T` by whatever type you use in the call". But if
944 there is no call later on in the same `.cc` file, then the compiler won't
945 do anything. If, at the same time, in a different `.cc` file that includes
946 the header file, you use the function with `T=double` the compiler will say
947 "Yes, I saw the declaration, but there is no definition; I assume the
948 function has been compiled in a different `.cc` file with `T=double` and
949 I'll simply record a call to this instantiation in the object file". The
950 call will then be resolved at link time if indeed another object file
951 contains an instantiation of the template for `T=double`. However, if no
952 other object file contains such a definition, a linker error will result.
953
954 In general, for functions like the above, it is difficult to foresee what
955 kinds of template arguments the function may be instantiated for, and so
956 there is no real practical way to put the definition into a `.cc` file.
957 Rather, one puts it into a header file, and so all `.cc` files that may use
958 this function see its definition (i.e. its body) and the compiler can
959 instantiate it in each source file for whatever template argument is
960 necessary. This makes sure that you never get linker errors, but at the
961 same time it makes compiling slow since every header file now not only has
962 to parse the function's declaration, but also its definition -- and in the
963 case of deal.II the definitions of all template functions add up to tens or
964 hundreds of thousands of lines of code. This is one of the reason why many
965 C++ programs compile relatively slowly: because they use a significant part
966 of the C++ standard library, most of which consists of templates.
967
968 deal.II can avoid much of this overhead. The trick is to recognize that in
969 the example above we don't really know what types `T` user code may
970 possibly want to use for this template. But in the case of using the space
971 dimension as a template parameter, we know pretty exactly all the possibly
972 values: `Triangulation<dim>` may really only be instantiated for `dim=1, 2,
973 3` and for nothing else. Consequently, we can do the following: Put all the
974 definitions of the member functions of deal.II into the `.cc` file and at
975 the bottom of the file instruct the compiler to please instantiate all of
976 these templates for `dim=1, 2, 3`. Similar things can be done for many
977 other template functions in deal.II; for example, there are a good number
978 of functions that require vector types as template arguments, of which
979 deal.II provides a good number, yet this list is finite and enumerable.
980 Consequently, we can simply, at the bottom of the `.cc` file, tell the
981 compiler to instantiate all of these template functions for every single
982 vector type deal.II supports, and then don't have to put thousands of lines
983 of template definitions into header files.
984
985 In many cases, enumerating all possible template arguments is tedious; it
986 is also difficult to extend this list when a new vector type is added, for
987 example. To simplify this task, deal.II uses a preprocessor: for many files
988 that want to instantiate a function or class for multiple template
989 arguments, we have a file `.inst.in` that has the equivalent of a
990 `for`-loop over all possible values or types for a template argument; the
991 file is processed by the `common/scripts/expand_instantiations` program to
992 produce a `.inst` file that can then be included into the `.cc` file.
993
994 ### Why do I need to use `typename` in all these templates?
995
996 This is indeed a frequent question. To answer it, it is necessary to
997 understand how a compiler deals with templates, which will take a bit of
998 space here. Let's take for example this case:
999
1000 ```cpp
1001   void f(int);
1002   void g(double d)
1003   {
1004     f(d);
1005   }
1006   void f(double);
1007 ```
1008
1009 Here, in the function `void g(double)`, we call `f` with a double as an
1010 argument. Because at that point the compiler has only seen the declaration
1011 of the first overload of `f`, it will convert the double `d` to an integer
1012 and call this first overload. The fact that a second overload was declared
1013 later does not change this situation, since it wasn't visible at the time
1014 the compiler parsed `g`.
1015
1016 Templates are designed to work essentially the same, but there are slight
1017 complications. Take this example:
1018
1019 ```cpp
1020   void f(int);
1021   void f(char);
1022   template <typename T> void g(T t)
1023   {
1024     f(1.1);
1025     f(t);
1026   }
1027   void f(double);
1028 ```
1029
1030 In the first line of `g`, the same thing happens as before: the argument is
1031 cast to `int` and the first of the two overloads of `f` is called. But when
1032 the compiler sees the template, it doesn't know yet what type `T` actually
1033 represents, so there is no way to settle on one of the two functions `f`
1034 the compiler has seen before when deciding about the second line. In fact,
1035 the C++ standard says that because the type of the argument `t` in the call
1036 depends on the template type, determining what function to actually call
1037 should only happen <i>at the time and place when the template is
1038 instantiated</i> (this is called <i>argument dependent name lookup</i> or
1039 <i>ADL</i>). In other words, if below the code above we had this:
1040 ```cpp
1041   void h()
1042   {
1043     g(1.1);
1044   }
1045 ```
1046
1047 then in the instantiation of `g` the first call would be to `f(int)`
1048 (because the argument 1.1 does not depend on the type given in the template
1049 argument, and consequently only functions are considered that were seen
1050 <i>before the definition of</i> `g(T)`) whereas the second call to `f`
1051 would be to `f(double)` -- even though `f(double)` wasn't even declared at
1052 the place the compiler saw the call in the template (though it is available
1053 at the place where we instantiate `g<double>`) -- because the function call
1054 argument `t` has type `T` and therefore depends on the template argument.
1055
1056 Argument dependent lookup allows you to use function templates like `g`
1057 with your own data types. For example, you could have your own library that
1058 does
1059 ```cpp
1060   #include <f.h>
1061
1062   struct X { /* something */ };
1063
1064   void f (const X & x) { /* do something with the X */ }
1065
1066   void my_function()
1067   {
1068     X x;
1069     g(x);
1070   }
1071 ```
1072
1073 Presumably the writer of the `g` function did not know about your own type
1074 `X` yet, but her code still works because you provided a suitable overload
1075 of `f` in your own code.
1076
1077 So ADL is clever and allows you to use templates in ways the author of the
1078 template did not anticipate. But it has a dark side: for every statement in
1079 your code, the compiler has to figure out whether it depends on the
1080 template types or not, and it needs in fact to know quite a lot about it.
1081 Take this example:
1082
1083 ```cpp
1084   int p;
1085   template <typename T>
1086   void g(T t)
1087   {
1088     T::something * p;
1089     f(p);
1090   }
1091 ```
1092
1093 Here, is the call to `f` dependent because `p` depends on the type `T`? If
1094 `f` is called with an argument of type `X` that is declared like this
1095
1096 ```cpp
1097   struct X
1098   {
1099     typedef int something;
1100   };
1101 ```
1102
1103 then `T::something * p;` would declare a local variable called `p` that is
1104 of type pointer-to-int. On the other hand, if we had
1105
1106 ```cpp
1107   struct X
1108   {
1109     static double something;
1110   };
1111 ```
1112
1113 then `T::something * p;` multiplies the variable `X::something` by the
1114 global variable `p` and ignores the result of the multiplication. The
1115 following call to `f` would then be non-dependent because the type of the
1116 (global) variable `p` does not depend on the template argument.
1117
1118 The example shows that the compiler can't know whether a call is dependent
1119 or not in a template it is just seeing unless we tell it that
1120 `T::something` is supposed to be a type or a variable or function name. To
1121 avoid this situation, C++ says: if a compiler sees `T::something` then this
1122 is a variable or function name unless it is prefixed by the keyword
1123 `typename` in which case it is supposed to be a type. In other words, the
1124 call to `f` here is going to be non-dependent:
1125 ```cpp
1126   int p;
1127   template <typename T>
1128   void g(T t)
1129   {
1130     T::something * p;
1131     f(p);
1132   }
1133 ```
1134
1135 and instantiating `g` with the first example for `X` is going to lead to
1136 errors because `T::something` didn't turn out to be a variable. On the
1137 other hand, if we had
1138 ```cpp
1139   int p;
1140   template <typename T>
1141   void g(T t)
1142   {
1143     typename T::something * p;
1144     f(p);
1145   }
1146 ```
1147
1148 then the call is dependent and will be deferred until the compiler knows
1149 the type of `T`.
1150
1151 ### Why do I need to use `this->` in all these templates?
1152
1153 This is a consequence of the same rule in the C++ standard as discussed in
1154 the previous question, Argument Dependent Lookup of names (ADL). Consider
1155 this piece of code:
1156 ```cpp
1157   template <typename T>
1158   class Base
1159   {
1160   public:
1161     void f();
1162   };
1163
1164   template <typename T>
1165   class Derived : public Base<T>
1166   {
1167   public:
1168     void g();
1169   };
1170
1171   template <typename T>
1172   void Derived<T>::g()
1173   {
1174     f();
1175   }
1176 ```
1177
1178 By the rules, when the compiler <i>parses</i> the function `Derived::g`
1179 (note that parsing happens before and independently of <i>instantiating</i>
1180 the function for a particular argument type `T`), it sees that the call to
1181 `f()` does not depend on the template type and so it looks for a
1182 declaration of such a function somewhere. In the example above, it doesn't
1183 find one (we'll come to this in a second), which will yield an error. On
1184 the other hand, in this code,
1185 ```cpp
1186   void f(); // global function
1187
1188   template <typename T>
1189   class Base
1190   {
1191   public:
1192     void f();
1193   };
1194
1195   template <typename T>
1196   class Derived : public Base<T>
1197   {
1198   public:
1199     void g();
1200   };
1201
1202   template <typename T>
1203   void Derived<T>::g()
1204   {
1205     f();
1206   }
1207 ```
1208
1209 it would find the global function and so when instantiating the function
1210 for, say, `T=int`, you'd get a function `Derived<int>::g` that would call
1211 the global function `::f`. This may or may not be what you had in mind.
1212
1213 The question of course is why the compiler didn't record a call to
1214 `Base<T>::f` in `Derived<int>::g`? After all, the compiler knows that
1215 `Derived` is derived from `Base`. This has a lot to do with the fact that
1216 at the time of <i>parsing</i> the template, the compiler doesn't know for
1217 which template arguments the template will later be instantiated, and with
1218 explicit or partial specializations.  Consider for example this code:
1219 ```cpp
1220   template <typename T>
1221   class Base
1222   {
1223   public:
1224     void f();
1225   };
1226
1227   class X
1228   {
1229   public:
1230     void f();
1231   };
1232
1233   template <> class Base<int> : public X {};
1234
1235   template <typename T>
1236   class Derived : public Base<T>
1237   {
1238   public:
1239     void g();
1240   };
1241
1242   template <typename T>
1243   void Derived<T>::g()
1244   {
1245     f();
1246   }
1247 ```
1248
1249 Here, if you look at `Derived<T>::g`, the call to `f()` will be resolved to
1250 `Base<T>::f` for all possible types `T`, unless `T=int` in which case the
1251 call will be to `X::f`. The point is that at the time the compiler sees
1252 (parses) the template, it simply doesn't know yet what `T` is, and so ADL
1253 says: if the call is not dependent, find a non-dependent function to record
1254 (e.g. a global function) rather than trying to find a call in scopes you
1255 can't yet know will be relevant (e.g. `Base` or `X`). Likewise, in this
1256 code,
1257 ```cpp
1258   template <typename T>
1259   class Base
1260   {
1261   public:
1262     void f();
1263   };
1264
1265   template <>
1266   class Base<int>
1267   {
1268   public:
1269     struct f {};
1270   };
1271
1272   template <typename T>
1273   class Derived : public Base<T>
1274   {
1275   public:
1276     void g();
1277   };
1278
1279   template <typename T>
1280   void Derived<T>::g()
1281   {
1282     f();
1283   }
1284 ```
1285
1286 the meaning of `f()` changes depending on the template type: if `T=int`, it
1287 creates an object of type `Base<int>>::f` and then throws the object away
1288 again immediately. For all other template arguments `T`, it calls
1289 `Base::f`.
1290
1291 Given this longish description of how compilers look up names under the ADL
1292 rule, let's get back to the original question: If you have this code,
1293 ```cpp
1294   template <typename T>
1295   class Base
1296   {
1297   public:
1298     void f();
1299   };
1300
1301   template <typename T>
1302   class Derived : public Base<T>
1303   {
1304   public:
1305     void g();
1306   };
1307
1308   template <typename T> void Derived<T>::g()
1309   {
1310     f();
1311   }
1312 ```
1313
1314 how do you achieve that the call in `Derived::g` goes to `Base::f`? The
1315 answer is: Tell the compiler to defer the decision of what the call is
1316 supposed to do till the time when it knows what `T` actually is. And we've
1317 already seen how to do that: we need to make the call <i>dependent</i> on
1318 `T`! The way to do that is this:
1319 ```cpp
1320   template <typename T>
1321   void Derived<T>::g()
1322   {
1323     this->f();
1324   }
1325 ```
1326
1327 Here, `this` is a pointer to an object of type `Derived<T>`, which is of
1328 course dependent. So the resolution of what the statement is supposed to
1329 represent is deferred until instantiation time; at that time, however, the
1330 compiler knows what the base class is (for example it knows if there are
1331 explicit specializations) and so it knows which base classes to look into
1332 in an attempt to find a function with the name `f`.
1333
1334 ### Does deal.II require C++11 support?
1335 The answer to this question depends on the version of deal.II that you are
1336 interested in using
1337
1338 #### deal.II version 9.0.0
1339 As of version 9.0.0, deal.II requires C++11 support equivalent to that provided
1340 by GCC 4.8.0., which is, essentially, every new feature in C++11.
1341
1342 #### deal.II version 8.5.0 and previous
1343 The current release of deal.II, 8.5.0, is compatible with the C++98 and C++03
1344 standards, but some features (e.g., the `LinearOperator` class) are only
1345 available if your compiler supports a subset of C++11 features. GCC 4.6 and
1346 newer implement enough of C++11 for these features to be turned on. More
1347 exactly, we currently require the following language features to be present:
1348
1349 1. `auto`-typed variables
1350 2. The `nullptr` keyword
1351 3. Move constructors
1352 4. The `declval` and `decltype` keywords
1353 5. Lambda functions
1354
1355 while we do not use the following features:
1356
1357 1. Marking virtual functions as `override`
1358 2. Some features of the `type_traits` header, such as `std::is_trivially_copyable`
1359 3. Inheriting constructors
1360 4. Template aliases
1361
1362 The deal.II documentation has a
1363 [page](http://dealii.org/8.5.0/doxygen/deal.II/group__CPP11.html) dedicated to
1364 the issue of what parts of C++11 we use and how this works. For a more complete
1365 list of features we do and do not use see
1366 [the GCC 4.6 C++11 compatibility page](https://gcc.gnu.org/gcc-4.6/cxx0x_status.html).
1367
1368 ### Can I convert Triangulation cell iterators to DoFHandler cell iterators?
1369
1370 Yes. You can also convert between iterators belonging to different
1371 DoFHandlers as long as the are based on the identical Triangulation:
1372
1373 ```cpp
1374 Triangulation<2>::active_cell_iterator it = triangulation.begin_active();
1375
1376 DoFHandler<2>::active_cell_iterator it2 (&triangulation, it->level(), it->index(), &dof_handler);
1377 ```
1378
1379 ## Questions about specific behavior of parts of deal.II
1380
1381 ### How do I create the mesh for my problem?
1382
1383 Before answering the immediate question, one remark: When you use adaptive
1384 mesh refinement, you definitely want the initial mesh to be as coarse as
1385 possible. The reason is that you can make it as fine as you want using
1386 adaptive refinement as long as you have memory and CPU time available.
1387 However, this requires that you don't waste mesh cells in parts of the
1388 domain where they don't pay off. As a consequence, you don't want to start
1389 with a mesh that is too fine to start with, because that takes up a good
1390 part of your cell budget already, and because you can't coarsen away cells
1391 that are in the initial mesh.
1392
1393 That said, there are essentially three ways to generate a mesh, all of
1394 which are discussed in significantly more detail in the step-49 tutorial
1395 program:
1396  - For many standard geometries (square, cube, circle, sphere, ...) there
1397    are functions in namespace `GridGenerator` that can generate coarse
1398    meshes.
1399  - If `GridGenerator` does not offer a mesh for the geometry you have, but
1400    if the geometry is simple, then you can often create one "by hand". Take
1401    a look, for example, at how we create the mesh in step-14 using the
1402    `Triangulation::create_triangulation` function. All you need to do is
1403    take a piece of paper, draw the geometry and a number of coarse cells
1404    that form quadrilaterals, identify the locations of vertices and the
1405    connectivity from cells to vertices, and pass the corresponding lists to
1406    the Triangulation. Something similar can be done for simple 3d
1407    geometries.
1408  - If your geometry is truly complicated enough so that you can't draw a
1409    mesh by hand any more (i.e. if it requires more than, for example, 20-30
1410    coarse mesh cells), then you'll need a mesh generator. For
1411    quadrilaterals and hexahedra, there aren't all that many mesh
1412    generators. [gmsh](http://www.gmsh.info),
1413    [lagrit](https://lagrit.lanl.gov/) and [cubit](http://cubit.sandia.gov/)
1414    come to mind. The primary problem is that most mesh generators' output
1415    meshes aren't particularly coarse by default, so you may want to pay
1416    particular attention to this point when running the mesh generator.
1417    (This is relevant since deal.II is particularly good about creating
1418    adaptively refined meshes, but if your coarse mesh is already very large
1419    then you will likely not have a lot of resources left to adaptively
1420    refine it some more.) Once you have a mesh from a mesh generator, you
1421    would read it using the `GridIn` class, as demonstrated, for example, in
1422    step-5.
1423  - As it was already mentioned, if you do need to work with a geometry for which all you have is a triangular or tetrahedral mesh, then you can convert this mesh into one that consists of quadrilaterals and hexahedra using the tethex program, see https://github.com/martemyev/tethex .
1424
1425 ### How do I describe complex boundaries?
1426
1427 You need to define classes derived from the `Boundary` base class and
1428 attach these to particular parts of the boundary of the triangulation. The
1429 `Triangulation` class will then query your boundary object whenever it
1430 needs a new point on the boundary after mesh refinement.
1431
1432 In deal.II releases after 8.1, the way geometry is described has been made
1433 much more flexible. In particular, it is no longer only possible to
1434 describe the boundary, but it is also possible to describe where points in
1435 the interior lie. The step-53 tutorial program explains how this is done
1436 for a realistic example. Additionally, the step-54 tutorial highlights how 
1437 boundary points can be mapped to CAD geometry.
1438
1439 #### How do I ensure that my refined grid is correct when using CAD geometry?
1440 When using CAD geometries, it is critical that the manifold ID's for objects 
1441 in the mesh (that is, faces and lines) are correctly set such that, upon refinement,
1442 the child elements are projected to the correct CAD entities. The mesh that
1443 results when this is not done correctly will not only be wrong, but may also
1444 difficult to interpret the source of the error.
1445
1446 The principles that constitute the current best practice are as follows:
1447 1. Start from the lowest codimension objects, and identify how to deform them. 
1448 For the case of a codimension one triangulation, the cells of a `Triangulation<2,3>` 
1449 are quadrilaterals, that should deform as a `TopoDS_FACE`. (For a pure 3d triangulation,
1450 you might start by setting the manifold IDs of the cells themselves.)
1451 2. For the codimension one case, attach your favourite manifold to the `TopoDS_Face` using 
1452 `cell->set_all_manifold_ids(manifold_id)`. Notice the use of `set_all_manifold_ids()`, 
1453 and **not** `set_manifold_id()`. This is necessary because you want all children 
1454 of these objects to belong to the same manifold. In particular, for this want all faces
1455 that are internal (i.e., that are shared between two objects with the same manifold 
1456 ID) to inherit the same manifold ID. For the choice of manifold to apply here,
1457 it is often most appropriate to employ a `NormalToMeshProjectionManifold`.
1458 3. Now go one codimension lower. For this codimension one example, that would mean setting
1459 the manifold IDs on curves (and for pure 3d meshes, surfaces). Follow the same rule as 
1460 before, setting all manifold IDs on faces that you know should follow a known 
1461 codimension one geometry (for curves, a known `TopoDS_EDGE`, or `TopoDS_WIRE`) using
1462 calls to `set_all_manifold_ids()`.  Here it is often most appropriate to choose 
1463 `ArcLengthProjectionManifold` objects, attached to the wires that identify your geometry. 
1464
1465 Doing things in this order guarantees that internal edges get the correct manifold ID.
1466
1467 ### How do I get the degree of freedom indices at vertices?
1468
1469 For simple cases (for example, where only `FE_Q` elements are used) you could
1470 use the `cell->vertex_dof_index()` function. This would mean that you'd require
1471 a single grid traversal to extract each DoF value  associated with support
1472 points corresponding to the vertices. The code to do this would look something
1473 like this:
1474 ```cpp
1475 auto cell = dof_handler.begin_active ();
1476 const auto endc = dof_handler.end ();
1477 for (; cell != endc; ++cell)
1478 {
1479   for (unsigned int vertex = 0; vertex < GeometryInfo<dim>::vertices_per_cell;
1480        ++vertex)
1481   {
1482     for (unsigned int component=0; component<fe.n_components(); ++component)
1483       // Index associated with the given component at the cell local vertex
1484       const unsigned int idx = cell->vertex_dof_index(vertex,component);
1485   }
1486 }
1487 ```
1488 If this is something that you might repeat often then you could use the `GridTools::find_cells_adjacent_to_vertex()` function to cache this association
1489 once up front for a single vertex, or the `GridTools::vertex_to_cell_map()` function
1490 to do the same for all vertices.
1491
1492 For the most general case (e.g., when using non-primitive finite elements), you
1493 might need to use `DoFTools::map_dofs_to_support_points()` function and then find
1494 which support points at located at the vertex position that you're interested in.
1495
1496 ### I am using discontinuous Lagrange elements (`FE_DGQ`) but they don't seem to have vertex degrees of freedom!?
1497
1498 Indeed. And here's the reason: a vertex is an entity that is shared between
1499 different cells, i.e. it doesn't belong to one cell or another. If you have
1500 a shape function that is associated with it, then its support will extend
1501 to all of the cells that are adjacent to the vertex since no cell is
1502 different than any other cell. This is what happens, for example, with the
1503 `FE_Q(1)` element. The same is true, by the way, for degrees of freedom
1504 (and associated shape functions) that correspond to edges and faces between
1505 cells.
1506
1507 But that doesn't answer the question of discontinuous elements. There, you
1508 have functions that are interpolation polynomials whose <i>support
1509 point</i> happens to be located at the same position as the vertex, but the
1510 actual support of the shape function is restricted to a single cell. In
1511 other words, '''logically''' these shape functions belong to a cell, not a
1512 vertex or edge or face, since the latter are all shared between adjacent
1513 cells. What this leads to is that, for example for the `FE_DGQ(1)` element,
1514 you have
1515  - `fe.dofs_per_vertex` is zero
1516  - `fe.dofs_per_line` is zero
1517  - `fe.dofs_per_face` is zero
1518  - `fe.dofs_per_cell` is 4 in 2d and 8 in 3d.
1519 In other words, all shape functions are associated with the cell interior.
1520
1521 If this answer isn't quite satisfactory (because, after all, the shape
1522 functions <i>are</i> defined by interpolation at the location of the
1523 vertices), one could turn the question around: If you ask me for the degree
1524 of freedom associated with vertex 13, then I should ask you in return
1525 <i>which one</i> you have in mind since if there, say, four cells that meet
1526 at this vertex, then there will be 4 degrees of freedom defined there.
1527 Likewise, if you ask me for the value of the degree of freedom associated
1528 with vertex 13, then I should ask you in return <i>which one</i> as the
1529 function is discontinuous there and will have multiple values at the
1530 location of the vertex.
1531
1532 ### How do I access values of discontinuous elements at vertices?
1533
1534 The previous question answered why DG elements aren't defined at the
1535 vertices of the mesh. Consequently, functions like `cell->vertex_dof_index`
1536 aren't going to provide anything useful. Nevertheless, there are occasions
1537 where one would like to recover values of a discontinuous field at the
1538 location of the vertices, for example to average the values one gets from
1539 all adjacent cells in recovery estimators.
1540
1541 So how does one do that? The answer is: Getting the values at the vertices
1542 of a cell works just like getting the values at any other point of a cell.
1543 You have to set up a quadrature formula that has quadrature points at the
1544 vertices and then use an FEValues object with it. If you then use
1545 FEValues::get_function_values, you will get the values at all quadrature
1546 points (i.e. vertices) at once.
1547
1548 Setting up this quadrature formula can be done in two different ways: (i)
1549 You can create an object of type `Quadrature` from a vector of points that
1550 you can initialize with the reference coordinates of the 2<sup>dim</sup>
1551 vertices of a cell; or (ii) you can use the `QTrapez` class that has its
1552 quadrature points in the vertices. In the latter case, however, you need to
1553 verify that the order of quadrature points is indeed the same order as the
1554 vertices of a cell and, if that is not the case, translate between the two
1555 numbering systems.
1556
1557 ### Does deal.II support anisotropic finite element shape functions?
1558
1559 There is currently no easy-to-use support for this. It's not going to work
1560 for continuous elements because we assume that `fe.dofs_per_face` is the
1561 same for all faces of a cell.
1562
1563 It may be possible to make this work for discontinuous elements, though.
1564 What you would have to do is define a bunch of different elements with
1565 anisotropic shape functions and select which element to use on which cells,
1566 using the `hp::DoFHandler` to deal with using different elements on
1567 different cells. The part that's missing is to implement elements with
1568 anisotropic shape functions. I imagine that this wouldn't be too
1569 complicated to do since the element is discontinuous, but someone would
1570 have to implement it.
1571
1572 That said, you can do anisotropic <i>refinement</i>, which of course also
1573 introduces a kind of anisotropic approximation of your finite element
1574 space.
1575
1576 ### The graphical output files don't make sense to me -- they seem to have too many degrees of freedom!
1577
1578 Let's assume you have a 2x2 mesh and a Q<sub>1</sub> element then you would
1579 assume that output files (e.g. in VTK format) just have 9 vertex locations
1580 and 9 values, one for each of the 9 nodes of the mesh. However, the file
1581 actually shows 16 vertices and 16 such values.
1582
1583 The reason is that frequently output quantities in deal.II are
1584 discontinuous: it may be that the finite element in use is discontinuous to
1585 begin with; or that the quantity we want to output is defined on a
1586 cell-by-cell basis (e.g. error indicators) and therefore discontinuous; or
1587 that it is a quantity computed from a DataPostprocessor object that could
1588 be discontinuous. In order to not make things more complicated than
1589 necessary, deal.II <i>always</i> assumes that quantities are discontinuous,
1590 even if some of them may in fact be continuous. The problem is that all
1591 graphical formats want to see one value for each output field per vertex.
1592 But discontinuous fields have more than one value at the location of a
1593 vertex of the mesh. The solution to the problem is then to simply output
1594 each vertex multiple times -- with different vertex numbers but at exactly
1595 the same location, once for each cell it is adjacent to. In other words, in
1596 2d, each cell has four unique vertices. The 2x2 mesh in the example
1597 therefore has 16 vertices (4 vertices for each of the 4 cells) and we
1598 output 16 values. Several of these vertices will have the same location and
1599 if the field is indeed continuous, several of the values will also be the
1600 same.
1601
1602 ### In my graphical output, the solution appears discontinuous at hanging nodes
1603
1604 Let me guess -- you are using higher order elements? If that's the
1605 case, then the solution only looks discontinuous but isn't
1606 really. What's happening is that the solution is, in fact, a higher
1607 order polynomial (e.g., a quadratic polynomial) along each edge of a
1608 cell but because all visualization file formats only support writing
1609 data as bilinear elements we need to write data in a way that shows
1610 only a linear interpolation of this higher order polynomial along each
1611 edge. This is no problem if the two neighboring elements share the
1612 entire edge because then the linear interpolations from both sides
1613 coincide. However, if we have a hanging node, then the value at the
1614 hanging node appears to float above or below the linear interpolation
1615 from the longer side, like here (in the left picture, see the gap at
1616 the bottom in the blue green area, and around the top left in the
1617 greenish area; pictures by Kevin Dugan):
1618
1619 <img width="400px" src="http://www.dealii.org/images/wiki/gap-in-q2-1.png" align="center" />
1620 <img width="400px" src="http://www.dealii.org/images/wiki/gap-in-q2-2.png" align="center" />
1621
1622 From this description you can already guess what the solution is: the
1623 solution is internally in fact continuous: even though we only show a
1624 linear interpolation on the long edge, the true solution actually goes
1625 through the "floating" node. All this is, consequently, just an
1626 artifact of the way visualization programs show data.
1627
1628 If this bothers you or it simply looks bad in your graphics, you can
1629 lessen the problem by not plotting just a linear interpolation on each
1630 cell but outputting the solution as a linear interpolation on a larger
1631 number of "patches" per cell (e.g., plotting 5x5 patches per
1632 cell). This can be done by using the `DataOut::build_patches` function
1633 with an argument larger than one -- see its documentation.
1634
1635 This all said, if you are in fact using a Q1 element and you see such
1636 gaps in the solution, then something is genuinely wrong. One
1637 possibility is that you forget to call `ConstraintMatrix::distribute`
1638 after solving the linear system, or you do not set up these
1639 constraints correctly. In either case, it's a bug if this happens with
1640 Q1 elements.
1641
1642 ### When I run the tutorial programs, I get slightly different results
1643
1644 This is sometimes unavoidable. deal.II uses a number of iterative
1645 algorithms (e.g. in solving linear systems, but the adaptive mesh
1646 refinement loop is also an iteration if you think about it) where certain
1647 criteria are specified by comparing floating point numbers. For example,
1648 the CG method terminates the iteration whenever the residual drops below a
1649 certain threshold; similarly, we refine as many cells as are necessary to
1650 take care of a fraction of the total error. In both cases, the quantities
1651 that are compared are floating point numbers which are subject to floating
1652 point round off. The problem is that floating point round off depends on
1653 the processor (sometimes), compiler flags or randomness (if parallelization
1654 is involved) and consequently an a solver may terminate one iteration
1655 earlier or later, depending on your environment, than the one from which we
1656 produced our results. With a different solution typically come different
1657 refinement indicators and different meshes downstream.
1658
1659 In other words, this is something that simply happens. What should worry
1660 you, however, is if you run the same program twice and you get slightly
1661 different output. This hints at non-deterministic effects that one should
1662 investigate.
1663
1664 ### How do I access the whole vector in a parallel MPI computation?
1665
1666 Note that this causes a bottleneck for large scale computations and you
1667 should try to use a parallel vector with ghost entries instead. If you
1668 really need to do this, create a TrilinosWrappers::Vector (or a
1669 PETScWrappers::Vector) and assign your parallel vector to it (or use a copy
1670 constructor). You can find this being done in step-17 if you search for
1671 "localized_solution".
1672
1673 ### How to get the (mapped) position of support points of my element?
1674
1675 Option 1: The support points on the unit cell can be accessed using
1676 FiniteElement::get_unit_support_point(s) and mapped to real coordinates
1677 using Mapping::transform_unit_to_real_cell()
1678
1679 Option 2: DoFTools::map_dofs_to_support_points() maps all the support
1680 points at once.
1681
1682 Option 3: You can create a FEValues object using the support points as a
1683 Quadrature:
1684 <pre>
1685 Quadrature<dim> q(fe.get_unit_support_points());
1686 FEValues<dim> fe_values (..., q, update_q_points);
1687 ...
1688 fe_values.get_quadrature_points();
1689 </pre>
1690
1691
1692 ## Debugging deal.II applications
1693
1694 ### I don't have a whole lot of experience programming large-scale software. Any recommendations?
1695
1696 Yes. First, the questions of this FAQ already give you a number of good
1697 pointers for example on debugging. Also, a good resource for some of the
1698 questions mathematicians, scientists and engineers (who may have taken a
1699 programming course, but know little of the bigger world of software
1700 engineering) typically have, is the [Software
1701 Carpentry](http://software-carpentry.org/) page. That site is specifically
1702 targeted at people who may want to use scientific computing to solve
1703 particular applications, but have little or no formal training in dealing
1704 with large software. In other words, it is specifically written for people
1705 for an audience like the users of deal.II.
1706
1707
1708 ### Are there strategies to avoid bugs in the first place?
1709
1710 Why yes, good you ask. There are indeed techniques that help you avoid
1711 writing code that has bugs. By and large, these techniques go by the name
1712 <i>defensive programming</i>, and the idea is to get yourself into a
1713 mindset while programming that anticipates that you will make mistakes,
1714 rather than expecting that your code is correct and then reacting to the
1715 situation when it turns out that this isn't true. The point is that even
1716 the most experienced programmers do introduce a lot of bugs into their
1717 code; what makes them good is that they have strategies to find them
1718 quickly and systematically.
1719
1720 Below we show one of the most important lessons learned. A more complete
1721 list can be found in [our code conventions
1722 page](http://dealii.org/developer/doxygen/deal.II/CodingConventions.html)
1723 which has a collection of best practices including code snippets to show
1724 how they are used.
1725
1726 The single most successful strategy to avoid bugs is to <i>make assumptions explicit</i>. For example, assume for a second that you have a class that denotes a point in 3d space:
1727 ```cpp
1728   class Point3d
1729   {
1730   public:
1731     double coordinate (const unsigned int i) const;
1732     // ...more here...
1733   private:
1734     double coordinates[3];
1735   };
1736
1737   double
1738   Point3d::coordinate (const unsigned int i) const
1739   {
1740     return coordinates[i];
1741   }
1742 ```
1743
1744 Here, when we wrote the `coordinate()` function, we worked under the
1745 assumption that the index `i` is between zero and two. As long as that
1746 assumption is satisfied, everything is fine. The problems start when
1747 someone calls this function with an index greater than two -- the function
1748 will in that case simply return garbage, but that may not be immediately
1749 obvious and may only much later lead to weird results in your program.
1750 Inexperienced programmers will say "Why would I do that, it doesn't make
1751 any sense!". Defensive programming starts from the premise that this is
1752 something that simply <i>will happen</i> at one point in time, whether you
1753 want to or not. It's actually not very difficult to do, since all of us
1754 have probably written code like this:
1755 ```cpp
1756   Point3d point;
1757   // ... do something with it
1758   double norm = 0;
1759   for (unsigned int i=0; i<=3; ++i)
1760     norm += point.coordinate(i) ** point.coordinate(i);
1761   norm = std::sqrt(norm);
1762 ```
1763
1764 Note that we have accidentally used `<=` instead of `<` in the loop.
1765
1766 If we accept that bugs will happen, we should make it as simple as possible
1767 to find them. In the spirit of making assumptions explicit, let's write
1768 above function like this:
1769 ```cpp
1770   double
1771   Point3d::coordinate (const unsigned int i) const
1772   {
1773     if (i >= 3)
1774       {
1775         std::cout << "Error: function called with invalid argument!" << std::endl;
1776         std::abort ();
1777       }
1778     return coordinates[i]
1779   }
1780 ```
1781
1782 This has the advantage that an error message is produced whenever the
1783 function is called with invalid arguments, and for good measure we also
1784 abort the program to make sure the error message can really not be missed
1785 in the rest of the output of the program. The disadvantage is that this
1786 check will always be performed whenever the program runs, even if it is
1787 well tested and we are fairly certain that in all places where the function
1788 is called, indices are valid. To avoid this drawback, the C programming
1789 language has the `assert` macro, which expands to the code above by
1790 default, but that can be disabled using a compiler flag. deal.II provides
1791 an improved version of this macro that is used as follows:
1792 ```cpp
1793   double
1794   Point3d::coordinate (const unsigned int i) const
1795   {
1796     Assert (i<3, ExcMessage ("Function called with invalid argument!"));
1797     return coordinates[i];
1798   }
1799 ```
1800
1801 The macro expands to nothing in optimized mode (see below), and if it is
1802 triggered in debug mode it doesn't only abort the program, but also prints
1803 an error message and shows how we got to this point in the program.
1804
1805 Using assertions in your program is the single most efficient way to make
1806 assumptions explicit and help find bugs in your program as early as
1807 possible. If you are looking for some more background, check out the
1808 wikipedia articles on
1809 [assertions](http://en.wikipedia.org/wiki/Assertion_(computing)),
1810 [preconditions](http://en.wikipedia.org/wiki/Precondition) and
1811 [postconditions](http://en.wikipedia.org/wiki/Postcondition), and generally
1812 the [design by contract
1813 methodology](http://en.wikipedia.org/wiki/Design_by_contract).
1814
1815 ### How can deal.II help me find bugs?
1816
1817 In addition to using the `Assert` macro introduced above, the deal.II
1818 libraries come in two flavors: debug mode and optimized mode. The
1819 difference is that the debug mode libraries contain a lot of assertions
1820 that verify the validity of parameters you may pass when calling library
1821 functions and classes; the optimized libraries don't contain these and are
1822 compiled with flags that instruct the compiler to optimize. This makes
1823 executables linked against the optimized libraries between 4 and 10 times
1824 faster. On the other hand, you will find that you will find 90% or more of
1825 your bugs by using the debug libraries because most bugs simply pass data
1826 to other functions that they don't expect or that don't make sense. The
1827 consequence is that you should always use debug mode when you are still
1828 developing your code. Only when it runs without bugs -- and under no
1829 circumstances any earlier -- should you switch to optimized mode to do
1830 production runs. One of the silliest things you can do is switch to
1831 optimized mode because you otherwise get an error you can't make sense of
1832 and that you don't know how to fix; certainly, if the library complains
1833 about something and you ignore it, nothing good can come out of the
1834 remainder of the run of your program.
1835
1836 You can switch between debug and optimized mode, at least for the example
1837 programs, by compiling the example with either `make debug` or `make
1838 release`. There are further
1839 [instructions](https://www.dealii.org/8.3.0/users/cmakelists.html#cmakesimple.build_type)
1840 in the documentation describing how to set this up in your own codes.
1841
1842 ### Should I use a debugger?
1843
1844 This question has an emphatic, unambiguous answer: Yes! You may get by for
1845 a while by just putting debug output into your program, compiling it, and
1846 running it, but ultimately finding bugs with a debugger is much faster,
1847 much more convenient, and more reliable because you don't have to recompile
1848 the program all the time and because you can inspect the values of
1849 variables and how they change. Learn how to use a debugger as soon as
1850 possible. It is time well invested.
1851
1852 Debuggers come in a variety of ways. On Linux and other Unix-like operating
1853 systems, they are almost all based in one way or other on the [GNU Debugger
1854 (GDB)](http://www.gnu.org/s/gdb/). GDB itself is a tool that is driven by
1855 interactively typing commands; if you know your way around with it, it is
1856 quite usable but it is rather austere and unless you are already familiar
1857 with this style of debugging, don't learn it. Rather, you should either use
1858 a graphical front-end or, even better, a front-end to GDB that is
1859 integrated into an Integrated Development Environment (IDE). An example of
1860 the stand-alone graphical front-ends to GDB are
1861 [DDD](http://www.gnu.org/software/ddd/), a program that was the first of
1862 its kind on Linux for many years but whose development has pretty much
1863 ceased in the early 2000s; it is still quite a good program, though.
1864 Another example is [KDbg](http://kdbg.org/), a GDB front-end for the KDE
1865 desktop environment.
1866
1867 As mentioned, a better choice is to use a debugger front-end that is
1868 integrated into the IDE. Every decent IDE has an integrated debugger, so
1869 you have your choice. A list of IDEs and how they work with deal.II is
1870 given in the C++ section of this FAQ.
1871
1872 ### deal.II aborts my program with an error message
1873
1874 You are likely seeing something like the following:
1875 ```
1876 --------------------------------------------------------
1877 An error occurred in line <1223> of file </.../dealii/include/deal.II/lac/vector.h> in function
1878     Number& dealii::Vector<Number>::operator()(dealii::Vector<Number>::size_type)
1879     [with Number = double; dealii::Vector<Number>::size_type = unsigned int]
1880 The violated condition was:
1881     i<vec_size
1882 The name and call sequence of the exception was:
1883     ExcIndexRangeType<size_type>(i,0,vec_size)
1884 Additional Information:
1885 Index 10 is not in the half-open range [0,10).
1886
1887 Stacktrace:
1888 -----------
1889 #0  ./deliberate-mistake: foo()
1890 #1  ./deliberate-mistake: main
1891 --------------------------------------------------------
1892
1893 ```
1894
1895 This error is generated by the following program:
1896 ```cpp
1897 #include <deal.II/lac/vector.h>
1898 using namespace dealii;
1899
1900 void foo ()
1901 {
1902   Vector<double> x(10);
1903   for (unsigned int i=0; i<=x.size(); ++i)
1904     x(i) = i;
1905 }
1906
1907
1908 int main ()
1909 {
1910   foo ();
1911 }
1912 ```
1913
1914 So what to do in a case like this? The first step is to carefully read what
1915 the error message actually says as it contains pretty much all the
1916 information you need. So let's take the error message apart:
1917
1918  - The first two lines tell you where the problem happened: in the current
1919    case, in line 1223 of file
1920    `/.../dealii/include/deal.II/lac/vector.h` in the function
1921    `Number& dealii::Vector<Number>::operator()(dealii::Vector<Number>::size_type)`.
1922    This is a function in the library, so you likely don't know what exactly it
1923    does and what to do with it, but there is more information to come.
1924
1925  - The second part is the condition that should have been true but wasn't,
1926    leading to the error: `i<vec_size`. The variables involved in this
1927    condition (`i,vec_size`) are local variables of the function, or member
1928    variables of the class, so again you may not be entirely familiar with
1929    them. But you can already gather some of the information: `i` likely is
1930    an index, which should have been less than the variable `vec_size`
1931    (which sounds a lot like the length of a vector); the assertion says
1932    that it <i>should</i> have been smaller, but that it wasn't actually.
1933
1934  - There is more information: The exception generated is of kind
1935    `ExcIndexRange<size_type>(i,0,vec_size)` and the additional information says
1936    `Index 10 is not in the half-open range [0,10)`. In other words, the variable
1937    `i` has value `10`, and `vec_size` is also ten. This should already give you
1938    a fairly good idea what is happening: the vector has size ten, and following
1939    C array convention, that means that only indices zero through nine are value,
1940    but ten is not.
1941
1942  - The final part of the error message -- the stack trace -- tells you how
1943    you got to this place: reading from the bottom, `main()` called `foo()`
1944    which called the function that generated the error.
1945
1946 Taken together, this information should allow you figure out in 80% of
1947 cases what was going on, and fix the problem. Here, it is that we used the
1948 condition `i<=x.size()` in the loop, rather than the correct condition
1949 `i<x.size()`. In the remaining 20% of cases, things might be more
1950 difficult. For example, `foo()` might be a large and difficult function,
1951 and you would need to know in which part of the function did we access an
1952 invalid index of the vector. Or `i` was an index computed from other
1953 variables and you'd need to find out why it got the invalid value. In these
1954 cases, you'll have to learn how to use a debugger such as gdb, and in
1955 particular how to move up and down in the call stack and to inspect local
1956 variables in your source code.
1957
1958 ### The program aborts saying that an exception was thrown, but I can't find out where
1959
1960 deal.II creates two kinds of exceptions (in deal.II language): ones where we
1961 simply abort the program because you are doing something that can't be right
1962 (such as accessing element 11 of a 10-element vector; this results in what has
1963 been discussed in the previous question) and ones that use the C++ construct
1964 `throw` to raise an exception. The latter construct is used for things that
1965 can't be statically checked in debug mode because they may depend on values
1966 read from input files or on a status that may simply change from one run of
1967 the program to the next; consequently, they <i>always</i> need to be verified,
1968 not only in debug mode, and there is sometimes a way to work around it in a
1969 program. The typical case is trying to write to a file that can't be opened
1970 (e.g. because the directory/file you specified in a parameter file doesn't
1971 exist or because the file system has run out of disk space).
1972
1973 Most of the time, the exceptions deal.II throws are annotated with the
1974 location and function where this exception was raised, and if you use a
1975 `main()` function such as the one used starting in step-6, this information
1976 will be printed. However, there are also cases where this kind of information
1977 is not available and then it is often difficult to establish where exactly the
1978 problem is coming from: all you know is that an exception was thrown, but not
1979 where or why.
1980
1981 To debug such problems, two approaches have proven useful:
1982
1983  - Run your program in a debugger (see the question about debuggers above,
1984    as well as these videos showing how to use the debugger in
1985    22c8e221823811aa1178b450171824af:
1986    http://www.math.tamu.edu/~bangerth/videos.676.8.html,
1987    http://www.math.tamu.edu/~bangerth/videos.676.25.html). You need to
1988    instruct the debugger to stop whenever an exception is thrown. If you
1989    work with gdb on the command line, then issue the command `catch throw`
1990    before starting the program and it will stop everytime the code executes
1991    a `throw` statement. Integrated development environments typically also
1992    have ways of switching this on. Note that not every exception that is
1993    thrown actually indicates an error -- sometimes, there are legitimate
1994    reasons to throw an exception and catch it in the calling function, so
1995    you may have to continue (resume) a number of times before finding the
1996    place where this happens.
1997
1998  - Debugging by subtraction: Starting at the end of your program, remove
1999    one function/code block after the other until your program runs through
2000    without aborting. For example, if your program looked like step-6, see
2001    if it runs through if you don't create graphical output in `run()`. If
2002    it does, then you know that the exception must have been thrown in the
2003    block of code you just removed. If the program continues to abort, then
2004    reduce the number of mesh refinement cycles to find out within which
2005    cycle the problem happens. If it happens in the very first cycle, then
2006    remove calling the linear solver. If the program now runs through, then
2007    the problem happened in the solver. If it still aborts, then it must
2008    have happened before the solver, for example in the assembly. Repeating
2009    this, you will be able to narrow down which statement caused the
2010    problem, and knowing where a problem happens is already more than half
2011    of what you need to know to fix it.
2012
2013
2014 ### I get an exception in `virtual dealii::Subscriptor::~Subscriptor()` that makes no sense to me!
2015
2016 The full text of the error message probably looks something like this (the
2017 stack trace at the bottom is of course different in your code):
2018 ```
2019 An error occurred in line <103> of file </.../deal.II/source/base/subscriptor.cc> in function
2020     virtual dealii::Subscriptor::~Subscriptor()
2021 The violated condition was:
2022     counter == 0
2023 The name and call sequence of the exception was:
2024     ExcInUse (counter, object_info->name(), infostring)
2025 Additional Information:
2026 Object of class N6dealii15SparsityPatternE is still used by 5 other objects.
2027   from Subscriber SparseMatrix
2028
2029 Stacktrace:
2030 -----------
2031 #0  /.../deal.II/lib/libdeal_II.g.so.7.0.0: dealii::Subscriptor::~Subscriptor()
2032 #1  /.../deal.II/lib/libdeal_II.g.so.7.0.0: dealii::SparsityPattern::~SparsityPattern()
2033 #2  /.../deal.II/lib/libdeal_II.g.so.7.0.0: dealii::BlockSparsityPatternBase<dealii::SparsityPattern>::reinit(unsigned int, unsigned int)
2034 #3  /.../deal.II/lib/libdeal_II.g.so.7.0.0: dealii::BlockSparsityPattern::reinit(unsigned int, unsigned int)
2035 #4  /.../deal.II/lib/libdeal_II.g.so.7.0.0: dealii::BlockSparsityPattern::copy_from(dealii::BlockCompressedSimpleSparsityPattern const&)
2036 #5  ./step-6: NavierStokesProjectionIB<2>::setup_system()
2037 #6  ./step-6: NavierStokesProjectionIB<2>::run(bool, unsigned int)
2038 #7  ./step-6: main
2039 ```
2040
2041 What is happening is this: deal.II derives a bunch of classes from the
2042 `Subscriptor` base class and then uses the `SmartPointer` class to point to
2043 such objects. `SmartPointer` is actually a fairly simple class: when given
2044 a pointer, it increases a counter in the `Subscriptor` base of the object
2045 pointed to by one, and when the pointer is reset to another object or goes
2046 out of scope, it decreases the counter again. (It can also records
2047 <i>who</i> points to this object.) If someone tries to delete the object
2048 pointed to, then the destructor `dealii::Subscriptor::~Subscriptor()` is
2049 run and checks that in fact the counter in this object is zero, i.e. that
2050 nobody is pointing to the object any more -- because if some pointer was
2051 still pointing to it, it would be a poor decision to delete the object as
2052 then the pointer would point to invalid memory. If the counter is nonzero,
2053 you get the error above: you are trying to delete an object that is still
2054 pointed to. In the case above, you try to delete a `SparsityPattern` object
2055 (that is, from the stack trace, a part of a block sparsity pattern) even
2056 though there is still a `SparseMatrix` pointing to it (we get this from the
2057 "Additional Information" field).
2058
2059 The solution in cases like these is to make sure that at the time you
2060 delete the object, no other objects still have pointers that point to it.
2061
2062 There is one rather frequent case that results in an error like the above
2063 and that is often difficult to understand: if an exception is thrown in
2064 some function and not caught, all local objects are destroyed in the
2065 opposite order of their declaration; if it isn't caught in the function
2066 that called the place where the exception was generated, its local
2067 variables are also destroyed, and so on. This automatic destruction of
2068 objects typically bypasses all the clean-up code you may have at the end of
2069 a function and can then lead to errors like the above. For example, take
2070 this code:
2071 ```cpp
2072 void f()
2073 {
2074   SparseMatrix s;
2075   SparsityPattern sp;
2076   // initialize sp somehow
2077   s.reinit (sp);
2078   Vector v;
2079   // build a linear system
2080
2081   solve_linear_system (s, v);
2082
2083   s.reinit ();
2084 }
2085 ```
2086
2087 If the code executes normally, at the bottom of the function, the local
2088 variables `s,sp,v` will be destroyed in reverse order. Since we have called
2089 `s.reinit()`, the object no longer stores a pointer to `sp` and so
2090 destruction of `sp` before `s` incurs no harm. But if the function
2091 `solve_linear_system` throws an exception, for example because the linear
2092 system is singular, the call to `s.reinit()` isn't executed any more, and
2093 you will get an error like the one shown at the top.
2094
2095 In cases like these, the challenge becomes finding where the exception was
2096 thrown. The easiest way is to run your program in a debugger and let the
2097 debugger tell you whenever an exception is generated. In `gdb`, you can do
2098 that by saying `catch throw` before running the program; essentially, the
2099 command puts a breakpoint on all places where exceptions are thrown.
2100 Remember, however, that not every place where an exception is thrown is a
2101 candidate for the problem above: it may also be an exception that is caught
2102 in the function above and that never propagates to a point where it
2103 produces trouble. Consequently, it may well happen that you have to
2104 continue several times after seeing an exception thrown until you finally
2105 find the place where the offending exception happens.
2106
2107 ### I get an error that the solver doesn't converge. But which solver?
2108
2109 Solvers are often deeply nested -- take a look for example at step-20 or
2110 step-22, where there is an outer solver for a Schur complement matrix, but
2111 both in the implementation of the Schur complement as well as in the
2112 implementation of the preconditioner we solve other linear problems which
2113 themselves may have to be preconditioned, etc. So if you get an exception
2114 that the solver didn't converge, which one is it?
2115
2116 The way to find out is to not wait till the exception propagates all the
2117 way to `main()` and display the error code there. Rather, you probably
2118 don't have a Plan B anyway if a solver fails, so you may want to abort the
2119 program if that happens. To do this, wrap the call to the solver in a
2120 try-catch block like this:
2121 ```cpp
2122   try
2123   {
2124     cg.solve (system_matrix, solution, system_rhs, preconditioner);
2125   }
2126   catch (...)
2127   {
2128     std::cerr << "*** Failure in Schur complement solver! ***" << std::endl;
2129     std::abort ();
2130   }
2131 ```
2132
2133 Of course, if this is in the Schur preconditioner, you may want to use a
2134 different error message. In any case, what this code does is catch the
2135 exceptions thrown by the solver here, or by the system matrix's `vmult`
2136 function (if not already caught there) or by the preconditioner (if not
2137 already caught there). If you had already caught exceptions in the `vmult`
2138 function and in the preconditioner, then you now know that any exception
2139 you get at this location must have been because the CG solver failed, not
2140 the preconditioner, etc. The upshot is that you need to wrap <i>every</i>
2141 call to a solver with such a try-catch block.
2142
2143 ### How do I know whether my finite element solution is correct? (Or: What is the "Method of Manufactured Solutions"?)
2144
2145 This is not always trivial, but there is an "industry-standard" way of
2146 verifying that your code works as intended, called the '''method of
2147 manufactured solutions'''. Before we describe the method, let us point this
2148 out: '''A code that has not been verified (i.e. for which correctness has
2149 not been established) is worthless. You do not want to have results in your
2150 thesis or a publication that may later turn out to be incorrect because
2151 your code does not converge to the correct solution!'''
2152
2153 The idea to verify a code is that you need a problem for which you know the
2154 exact solution. Unless you solve the very simplest possible partial
2155 differential equations, it is typically not possible to choose a right hand
2156 side and boundary values and then find the corresponding solution to the
2157 PDE analytically, on a piece of paper. But you can turn this around: Let's
2158 say your equation is <i>Lu=f</i>, then choose some function <i>u</i> and
2159 compute <i>f=Lu</i>. Note that the solution <i>u</i> does not necessarily
2160 have to be something that looks like a useful or physically reasonable
2161 solution to the equation, all that is necessary is that it is a function
2162 you know. Because <i>L</i> is a differential operator, computing <i>f</i>
2163 only involves computing the derivatives of the known function; this may
2164 yield lengthy expressions if you have nonlinearities or spatially variable
2165 coefficients in the equation, but should not be too complicated and can
2166 also be done using computer algebra programs such as Maple or Mathematica.
2167
2168 If you now put this particular right hand side <i>f</i> into your program
2169 (along with boundary values that correspond to the values of the function
2170 <i>u</i> you have chosen) you will get a numerical solution
2171 <i>u<sub>h</sub></i> that we would hope converges against the exact
2172 solution <i>u</i> at a particular rate, say <i>O(h<sup>2</sup>)</i> in the
2173 <i>L<sub>2</sub></i> norm. But since you know the exact solution (you have
2174 chosen it before), you can compute the error between numerical solution and
2175 exact solution, and verify not only that your code converges, but also that
2176 it shows the convergence rate you expect.
2177
2178 The method of manufactured solutions is shown in the step-7 tutorial
2179 program.
2180
2181 ### My program doesn't produce the expected output!
2182
2183 There are of course many possible causes for this, and you need to find out
2184 which of these causes might be the reason. Possible places to start are:
2185
2186  - Are matrix and right hand side assembled correctly? For most reasonably
2187    simple problems, you can compute the local contributions to these
2188    matrices by hand, and then compare those with the ones you compute on
2189    every cell of your program (remember that you can print the contents of
2190    the local matrix and right hand side to screen). A good strategy is also
2191    to reduce your problem to a 1x1 or 2x2 mesh and then print out the
2192    entire system matrix for comparison.
2193
2194  - Do you compute the matrix you need, or its transpose? The mathematical
2195    literature often multiplies the equation from the right with a test
2196    function but that is awkward because the matrix you get this way is the
2197    transpose from the one you need. The deal.II documentation goes to
2198    lengths in multiplying test functions from the left to avoid this sort
2199    of error; do the same in your derivations.
2200
2201  - Your constraints or boundary values may be wrong. While the
2202    ConstraintMatrix and functions like
2203    VectorTools::interpolate_boundary_values are well enough tested that
2204    they are unlikely candidates for problems, you may have computed
2205    constraints wrongly if you collect them by hand (for example if you deal
2206    with periodic boundary conditions or similar) or you may have specified
2207    the wrong boundary indicator for a Dirichlet boundary condition. Again,
2208    the solution is to reduce the problem to the simplest one you can find
2209    (e.g. on the 1x1 or 2x2 mesh talked about above) and to ask the
2210    ConstraintMatrix to print its contents so that you can compare it by
2211    hand with your expectations.
2212
2213  - Your discretization might be wrong. Some equations require you to use
2214    particular (combination of) finite elements; for example, for the Stokes
2215    equations and many other saddle point problems, you need to satisfy an
2216    LBB or Babuska-Brezzi condition. For other equations, you need to add
2217    stabilization terms to the bilinear form; for example, advection or
2218    transport dominated problems require stabilization terms such as
2219    artificial diffusion, streamlinear diffusion, or SUPG.
2220
2221  - The solver might be wrong. This can reasonably easily happen if you have
2222    a complex solver such as, for example, the one used in step-22. In such
2223    cases it has proven useful to simply replace the entire solver by the
2224    sparse direct UMFPACK solver (see step-29). UMFPACK is not the fastest
2225    solver around, but it never fails: if the linear system has a solution,
2226    UMFPACK will find it. If the output of your program is essentially the
2227    same as before, then the solver wasn't your problem.
2228
2229  - Your assumptions may be wrong. Double check that you had the correct
2230    right hand side to compute the numerical solution you compare against
2231    your analytical one. Also remember that the numerical solution is
2232    usually only an approximation of the true one.
2233
2234 In general, if your program is not computing the output you expect, here
2235 are a few strategies that have often worked in finding the problem:
2236  - Take a good look at the output you get. For example, a close look can
2237    already tell you if (i) the boundary conditions are correct, (ii) the
2238    solution is continuous at hanging nodes, (iii) the solution follows the
2239    characteristics of the right hand side. This may already help you narrow
2240    down which part of the program may be the culprit. A common mistake is
2241    also to have a solution that by some accident is too large by a certain
2242    factor; consequently, the error will not converge to zero but to some
2243    constant value. This, again, is easily visible from a graphical
2244    representation of the solution and/or the error. Plotting the error is
2245    discussed in the section below entitled "How to plot the error as a
2246    pointwise function".
2247  - If you have a time dependent problem, is the first time step right?
2248    There is no point in running the program for 1000 time steps and trying
2249    to find our why it is wrong, if already the first time step is wrong.
2250  - If you still can't find what's going on, make the program as small as
2251    possible. Copy it to another directory and start stripping off parts
2252    that you don't need. For example, if it is a time dependent program for
2253    which you have previously already found out that the first time step is
2254    wrong, then remove the time loop. If you have tried whether you have the
2255    same problem when the mesh is uniformly refined, then throw out all the
2256    code that deals with adaptive refinement, constraints and hanging nodes.
2257    In this process, every time you simplify the program, verify that the
2258    problem is still there. If the problem disappears, you know that it must
2259    have been in the last simplification step. If the problem remains, it
2260    must be in the code that is now one step smaller. Ultimately, the code
2261    should be small enough so that you can just go through it and find the
2262    error by inspection.
2263  - Learn to use a debugger. You will find that using a debugger is so much
2264    more convenient than trying to put screen output statements into your
2265    code, recompiling, and hoping that they reveal the problem. Modern
2266    integrated development environment as the ones discussed elsewhere in
2267    this FAQ have the debugger built-in, allowing you to use it seamlessly
2268    in your editing environment.
2269
2270 ### The solution converges initially, but the error doesn't go down below 10<sup>-8</sup>!
2271
2272 First: If the error converges to zero, then you are basically doing
2273 something right already. Congratulations!
2274
2275 As for why the error does not converge any further, there are two typical
2276 cases what could be the reason:
2277
2278  - While the discretization error should converge to zero, the error of
2279    your numerical solution is composed of both the discretization error and
2280    the error of your linear or nonlinear solver. If, for example, you solve
2281    the linear system to an accuracy of 10<sup>-5</sup>, then there will be
2282    a point where the discretization error will get smaller than that by
2283    using finer and finer meshes but the solver error will not become
2284    smaller any more. To continue observing the correct convergence order,
2285    you will also have to solve the linear system with more accuracy.
2286
2287  - If you compute the error through an external program, for example by
2288    writing out the solution to a file and reading it from another program
2289    that knows about the exact solution, then you need to make sure you
2290    write the solution with sufficient accuracy. The default setting of C++
2291    writes floating point numbers with approximately 8 digits, so if you
2292    want to make sure that your solution is correct to 10<sup>-10</sup>, for
2293    example, you'll have to write out the solution with more than 10 digits.
2294
2295
2296 ### My code converges with one version of deal.II but not with another
2297
2298 That is a tough case because the problem could literally be anywhere in the functions you call from deal.II.
2299 Rather than trying to start debugging blindly to find out what exactly is going on it's probably more productive to delineate the steps one could use to narrow down where the problem is.
2300
2301 In an ideal world, you would have already found out which commit in the history of deal.II caused the problem.
2302 Let's say you have checked out the two offending versions of deal.II into separate source directories `dealii-good` and `dealii-bad`, and that you compiled them both separately and installed them into directories `install-good` and `install-bad`. If you can't find out which commit caused the problem, the good and bad versions could also be the last two releases.
2303
2304 Let's also say that you have a directory `application` in which you have your own code.
2305 Now create two directories, `app-good`, `app-bad` parallel to `application`. Then do
2306 ```bash
2307   cd application
2308   for i in * ; do
2309     ln -s $i ../app-good/$i
2310     ln -s $i ../app-bad/$i
2311   done
2312 ```
2313 This way you have two directories in which you can configure, compile, and run the exact same version of your application (exact same because they both contain links to the exact same source files), just compiled against the good and bad versions of the library, respectively.
2314
2315 So you do
2316 ```bash
2317   cd app-good
2318   cmake . -DDEAL_II_DIR=.../install-good
2319   make
2320
2321   cd ../app-bad
2322   cmake . -DDEAL_II_DIR=.../install-bad
2323   make
2324 ```
2325 If you run in these two directories, e.g., in two separate xterm windows, you will get one working and one failing run. Now start modifying the source files in `application` to figure out where the first point in the program is where there are differences. For example, after assembly, you could do insert a statement of the form
2326 ```cpp
2327   std::cout << "Linear system: " << system_matrix.l1_norm() << ' ' << system_rhs.l2_norm() << std::endl;
2328 ```
2329 I would suspect (though that doesn't have to be true -- but just assume for the moment) that if you compile and run this modification in your two windows that you will get different results. At this point, you can remove everything that is executed after this point from your program -- likely a few hundred lines of code. Or, if you're too lazy, just put `abort()` after that statement because everything that comes after it clearly only shows symptoms but not the cause of the problem.
2330
2331 Now that you know that the problem exists at the end of assembly, make your way further forward in the program. For example, is the local matrix on the first cell on which you assemble the same between the two programs? If it is, the problem is on a later cell. If it isn't the same, try to think about what the cause may be. Is the mesh the same? You can test that by putting output into an earlier spot of the program; if that output is different between the two programs, you can again delete everything that happens after that point.
2332
2333 The whole exercise is designed to find the first place in the program where you can unambiguously say that something has changed. Non-convergence is just such a non-specific problem that it is not helpful in finding what exactly is going on. It also happens rather late in typical programs that there are too many possibilities for where the root cause may be.
2334
2335 ### My time dependent solver does not produce the correct answer!
2336
2337 For time dependent problems, there are a number of other things you can try
2338 over the discussion already given in the previous answer. In particular:
2339
2340  - If you have a time iteration and the solution at the final time (where
2341    you may evaluate the error) is wrong, then it was likely already wrong
2342    at the first time step. Try to run your program only for a single time
2343    step and make sure the solution there is correct. For example, it could
2344    be that you set the boundary values wrongly; this would be quite
2345    apparent if you looked at the first time step because the effect would
2346    be largest close to the boundary, but it may no longer be visible if you
2347    ran your program for a couple hundred time steps.
2348
2349  - Are your initial values correct? Output the initial values using DataOut
2350    just like you output the solution and inspect it for correctness.
2351
2352  - If you have a multi-stage time stepping scheme, are *all* the initial
2353    values correct?
2354
2355  - Finally, you can test your scheme by setting the time step to zero. In
2356    that case, the solution at time step zero should of course be equal to
2357    the solution at time step zero. If it isn't, you already know better
2358    where to look.
2359
2360 ### My Newton method for a nonlinear problem does not converge (or converges too slowly)!
2361
2362 Newton methods are tricky to get right. In particular, they sometimes
2363 converge (if slowly) even though the implementation has a bug because all
2364 that is required for convergence is that the search direction is a
2365 direction of descent; consequently, if for example you have the wrong
2366 matrix, you may compute something that is a direction of descent, but not
2367 the full Newton direction, and so converges but not at quadratic order.
2368
2369 Here are a few considerations for implementing Newton's method for
2370 nonlinear PDEs:
2371
2372  - Try it with a linear program by removing all the nonlinearities in your
2373    problem. Your Newton iteration must converge in a single step, i.e. the
2374    Newton residual must be zero at the beginning of the second iteration.
2375    If that's not the case, something is wrong in your implementation.
2376
2377  - Newton's iteration will converge with optimal order for the problem
2378    *R(u)=0*, where *R* may be thought of as a residual, if you
2379    <i>consistently</i> compute the Newton residual
2380    *(&phi;<sup>i</sup>, R(u<sup>k</sup>))* and the Newton (Jacobian) matrix
2381    *R'(u<sup>k</sup>)*. If you have
2382    a bug in either of the two, your method may converge, but typically at a
2383    (much) lower rate and with consequently many more iterations.
2384
2385    Consequently, one way to debug Newton's methods is to verify that the Newton
2386    matrix and Newton residual are matching in their code. However, if you
2387    have a matching bug in <i>both</i> of the matrix and right hand side
2388    assembly, then your Newton method will converge with correct order but
2389    against the wrong solution.
2390
2391  - If you have nonzero boundary values for your problem, set the correct
2392    boundary values for the initial guess and use zero boundary values for
2393    all following updates. This way, the updated
2394    *u<sup>k+1</sup> = u<sup>k</sup> + &delta; u<sup>k</sup>*
2395    already has the right boundary values for all following iterations, where
2396    *&delta; u<sup>k</sup>* is the Newton update.
2397
2398  - If your problem is strongly nonlinear, you may need to employ a line
2399    search where you compute
2400    *u<sup>k+1</sup> = u<sup>k</sup> + &alpha; &delta; u<sup>k</sup>*
2401    and successively try *&alpha;=1, &alpha;=1/2, &alpha;=1/4*, etc., until the
2402    residual computed for *u<sup>k+1</sup>* for this *&alpha;* is smaller than
2403    the residual for *u<sup>k</sup>*.
2404
2405  - A rule of thumb is that if your problem is strongly nonlinear, you may
2406    need 5 or 10 iterations with a step length *&alpha;* less than one, and
2407    all following steps use the full step length *&alpha;=1*.
2408
2409  - For most reasonably behaved problems, once your iteration reaches the
2410    point where it takes full steps, it usually converges in 5 or 10 more
2411    iterations to very high accuracy. If you need significantly more than 10
2412    iterations, something is likely wrong.
2413
2414 ### Printing deal.II data types in debuggers is barely readable!
2415
2416 Indeed. For example, plain gdb prints this for cell iterators:
2417 ```
2418 $2 = {<dealii::TriaIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 3> > >> = {<dealii::TriaRawIterator<dealii::DoFCellAccessor<dealii::DoFHandler<2, 3> > >> = {<std::iterator<std::bidirectional_iterator_tag, dealii::DoFCellAccessor<dealii::DoFHandler<2, 3> >, long, dealii::DoFCellAccessor<dealii::DoFHandler<2, 3> >*, dealii::DoFCellAccessor<dealii::DoFHandler<2, 3> >&>> = {<No data fields>},
2419       accessor = {<dealii::DoFAccessor<2, dealii::DoFHandler<2, 3> >> = {<dealii::CellAccessor<2, 3>> = {<dealii::TriaAccessor<2, 2, 3>> = {<dealii::TriaAccessorBase<2, 2, 3>> = {
2420                 static space_dimension = <optimized out>, static dimension = <optimized out>,
2421                 static structure_dimension = <optimized out>, present_level = -9856,
2422                 present_index = 32767, tria = 0x4a1556}, <No data fields>}, <No data fields>},
2423           static dimension = 2, static space_dimension = 3, dof_handler = 0x7fffffffdac8},
2424         static dim = <optimized out>,
2425         static spacedim = <optimized out>}}, <No data fields>}, <No data fields>}
2426 ```
2427
2428 Fortunately, this can be simplified to this:
2429 ```cpp
2430 $3 = {
2431   triangulation = 0x4a1556,
2432   dof_handler = 0x7fffffffdac8,
2433   level = 2,
2434   index = 52
2435 }
2436 ```
2437
2438 All you need is (i) gdb version 7.1 or later, or a graphical frontend for
2439 it (e.g. [DDD](http://www.gnu.org/software/ddd/) or
2440 [kdevelop](http://www.kdevelop.org)), (ii) some code that goes into your
2441 $HOME/.gdbinit file. Instructions for setting up this file, which implements
2442 pretty printers for `Point`, `Tensor`, `Vector`, and the various iterator
2443 classes for triangulations and DoFHandlers, are posted
2444 [here](Debugging-with-GDB).
2445
2446 gdb can also pretty print many of the `std::XXX` classes, but not all linux
2447 distributions have it configured this way. To enable this, follow the
2448 instructions [from this
2449 website](http://sourceware.org/gdb/wiki/STLSupport). The little python
2450 snippet can be placed as a separate python block into `.gdbinit`.
2451
2452 ### My program is slow!
2453
2454 This is a problem that is true for a lot of us. The question is which part
2455 of your program is causing it. Before going into more detail, there are,
2456 however, some general observations:
2457
2458  - Running deal.II programs in debug mode will take, depending on the program,
2459    between 4 and 10 times as long as in optimized mode. If you are using the
2460    standard setup for your own `CMakeLists.txt` file (described in the
2461    [documentation](https://www.dealii.org/8.3.0/users/cmakelists.html)), then
2462    compiling your code with `make release` will both compile your code at a
2463    higher optimization level and link it against the optimized version of
2464    deal.II.
2465
2466  - A typical finite element program will spend around one third of its time
2467    in assembling linear systems, around one half in solving these linear
2468    systems, and the rest of the time on other things. If your program's
2469    percentages significant deviate from this rule of thumb, you know where
2470    to start.
2471
2472  - There is a rule that says that even the best programmers are unable to
2473    point out where in the program the most CPU time is spent without some
2474    form of profiling. This is definitely true also for the primary
2475    developers of deal.II, so it is likely true for you as well. A corollary
2476    to this rule is that if you start optimizing parts of your code without
2477    first profiling it, you are more than likely just going to make things
2478    more complicated without significant gains because you pick the simplest
2479    places to optimize, not the ones with the biggest impact.
2480
2481 So how can you find out which parts of the program are slow? There are two
2482 tools that we've really come to like, both from the
2483 [valgrind](http://www.valgrind.org/) project: callgrind and cachegrind.
2484 Valgrind essentially emulates what your CPU would do with your program and
2485 in the process collects all sorts of information. In particular, if you run
2486 your program as in
2487 ```
2488   valgrind --tool=callgrind ./myprogram
2489 ```
2490
2491 (this will take around 10 times longer than when you just call
2492 `./myprogram` because of the emulation) then the result will be a file in
2493 this directory that contains information about where your program spent its
2494 time. There are a number of graphical frontends that can visualize this
2495 data; my favorite is `kcachegrind` (a misnomer -- it is, despite its name,
2496 actually a frontend from callgrind, not cachegrind). Pictures of how this
2497 output looks can be found in the introduction of step-22. It typically
2498 shows how much time was spent in each function and a call graph of which
2499 functions where called from where.
2500
2501 Using valgrind's cachegrind can give you a more detailed look at much of
2502 the same kind of information. In particular, it can show you source line
2503 for source line how many instructions were executed there, and how many
2504 memory accesses (as well as cache hits and misses) were generated there.
2505 See the valgrind manual for more information.
2506
2507 Lastly, since you are probably most interested in the performance of the
2508 optimized version of your code (which you will probably use for long
2509 expensive runs), you should run valgrind on the optimized executable.
2510
2511 ### How do I debug MPI programs?
2512
2513 This is clearly an awkward topic for which there are few good options:
2514 debugging parallel programs using MPI has always been a pain and it is
2515 frustrating even to experienced programmers. That said, there are parallel
2516 debuggers that can deal with MPI, for example
2517 [TotalView](http://www.roguewave.com/products/totalview-family/totalview.aspx)
2518 that can make this process at least somewhat simpler.
2519
2520 Whether you have or don't have TotalView, here are a few guidelines of
2521 strategies that have helped us in the past:
2522
2523  - Try to reduce the problem to the smallest one you can find: The smallest
2524    mesh, the smallest number of processors. Reducing the number of
2525    processors needed to demonstrate the bug must be your highest priority.
2526
2527  - One of the biggest problems you typically have is that the processes
2528    that communicate via MPI typically run on different machines. If you can
2529    manage to reduce the problem to a small enough number of processors, you
2530    can run them all locally on a single workstation, rather than a cluster
2531    of computers. Ideally, you would reduce the problem to 2 or 4 processors
2532    and then just start the program using `mpirun -np 4 ./myexecutable` on
2533    the headnode of the cluster, a workstation, or even a laptop.
2534
2535  - Try to figure out which MPI process (the MPI rank) has the problem, for
2536    example by printing the output of
2537    `Utilities::System::get_this_mpi_process(MPI_COMM_WORLD)` at various
2538    points in your program.
2539
2540  - If you know which MPI process has the problem and if this is
2541    reproducible, let each process print out its MPI rank and its process id
2542    (PID) using the system function `getpid` at the very beginning of the
2543    program. The PID is going to be different every time you run the
2544    program, but if you know the connection between MPI rank and PID and you
2545    know which rank will produce the problem, then you can predict which PID
2546    will have the problem. The point of this is that you can attach a
2547    debugger to this PID; for example, `gdb` has the command `attach <pid>`
2548    with which you can attach the debugger to a running program, rather than
2549    running the program from the start within the debugger. Attaching a
2550    program to the debugger will stop it (and, after a while, will typically
2551    also stop all the other MPI processes once they come to a place where
2552    they are waiting for a communication from the stopped process). You can
2553    then look at variables, continue running to breakpoints, or do whatever
2554    else you want to do with the process you attached the debugger to. In
2555    particular, if for example you attached the debugger to the process that
2556    you know will segfault or run onto a failing assertion, you can just
2557    type `continue` in the debugger to let the program continue till it
2558    aborts. You can then inspect the state of the program at the point of
2559    the problem inside the debugger you attached.
2560
2561  - The above process relies on the fact that you have time to attach a
2562    debugger between starting the program, reading the mapping from MPI
2563    process rank to PID, and attaching a debugger. If the program produces
2564    the error very quickly, it is often useful to insert a call to
2565    `sleep(60);` (and including the appropriate header file) just after
2566    outputting MPI rank and PID. This gives you 60 seconds to attach the
2567    debugger before the program will continue.
2568
2569  - If finding out which MPI process has the problem turns out to be too
2570    complicated, or if it isn't predictable which process will produce an
2571    error, then there is a fallback option: attach a debugger to
2572    <i>every</i> MPI process. This is awkward to do by hand, but there is a
2573    shortcut: under linux (or any other unix system) you can run
2574    the program as in
2575 ```
2576   mpirun -np 4 xterm -e gdb --args ./my_executable
2577 ```
2578 The equivalent for macOS is
2579 ```
2580   mpirun -np 4 xterm -e lldb -f ./my_executable
2581 ```
2582
2583 In this example, we start 4 MPI processes; in each of these 4 processes, we
2584 open an `xterm` window in which we start an instance of `gdb`/`lldb` with the
2585 executable. You'd then `run` the executable in each of the 4 windows, and
2586 debug it as you usually would. This might be tedious but as mentioned
2587 above, debugging MPI programs often is tedious indeed. To find out which
2588 gdb window belongs to which MPI rank, you can type the command
2589 ```
2590   !env | grep RANK
2591 ```
2592 into the gdb window (this works with OpenMPI at least).
2593
2594 ### I have an MPI program that hangs
2595
2596 Apart from programs that segfault or that run onto a failing assertion
2597 (both cases that are relatively easy to debug using the techniques above),
2598 programs that just hang are the most common problem in parallel
2599 programming. The typical cause for this is that there is a point in your
2600 program where all or some MPI processes expect to get a message from a
2601 process X (e.g. in a global communication, say MPI_Reduce, MPI_Barrier, or
2602 directly in point-to-point communications) but process X is not where it
2603 should be -- for example, because it is in an endless loop, or -- more
2604 likely -- because process X didn't think that it should participate in this
2605 communication. In either case, the other processes will wait forever for
2606 process X's message and deadlock the program. An example for this case
2607 would go like this:
2608 ```cpp
2609   void assemble_system ()
2610   {
2611     // optimization in case there is nothing to do; we won't
2612     // have to initialize FEValues and other local objects in
2613     // that case
2614     if (tria.n_locally_owned_active_cells() == 0)
2615       return;
2616
2617     ...
2618     for (cell = ...)
2619       if (!cell->is_ghost() && !cell->is_artificial())
2620         ...do the assembly on the locally owned cells...
2621
2622     system_rhs.compress();
2623   }
2624 ```
2625
2626 Here, the call to `compress()` at the end involves communication between
2627 MPI processes. In particular, say, it implies that process Y will wait for
2628 some data from process X. Now what happens if process X realizes that it
2629 doesn't have any locally owned cells? In that case, process X will quit the
2630 function at the very top, and will never call `compress()`. In other words,
2631 process Y will wait forever, possibly making process Z wait further down
2632 the program etc. In the end, the program will be deadlocked.
2633
2634 The goal of debugging the program must be to find where individual
2635 processes are stopped in order to determine which incoming communication
2636 they are waiting for. If you attached a debugger to the program above,
2637 you'd find for example that all but one process is stopped in the call to
2638 `compress()`, and the one remaining process is stopped in some other MPI
2639 call, then you already have a good idea what may be going on.
2640
2641 ### One statement/block/function in my MPI program takes a long time
2642
2643 Let's say you have a block of code that you suspect takes a long time and
2644 you want to time it like this:
2645 ```cpp
2646   Timer t;
2647   t.start();
2648   my_function();
2649   t.stop();
2650   if (my MPI rank == 0)
2651     std::cout << "Calling my_function() took " << timer() << " seconds." << std::endl;
2652 ```
2653
2654 The output is large, i.e. you think that the function you called is taking
2655 a long time to execute and that you should focus your efforts on optimizing
2656 it. But in an MPI program, this isn't quite always true. Imagine, for
2657 example, that the function looked like this:
2658 ```cpp
2659   void my_function ()
2660   {
2661     double val = compute_something_locally();
2662     double global_sum = 0;
2663     MPI_Reduce (&val, &global_sum, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
2664     if (my MPI rank == 0)
2665       std::cout << "Global sum = " << global_sum << std::endl;
2666   }
2667 ```
2668
2669 In the call to `MPI_Reduce`, all processors have to send something to
2670 processor zero. Processor zero will have to wait till everyone sends stuff
2671 to this processor. But what if processor X is still busy doing something
2672 else (stuff above the call to `my_function`) for a while? The processor
2673 zero will wait for quite a while, not because the operations in
2674 `my_function` are particularly expensive (either on processor zero or
2675 processor X) but because processor X was still busy doing something else.
2676 In other words: you need to direct your efforts in making the "something
2677 else on processor X" faster, not making `my_function` faster.
2678
2679 To find out whether this is really the problem, here is a simple way to see
2680 what the "real" cost of `my_function` is:
2681 ```cpp
2682   Timer t;
2683   MPI_Barrier (MPI_COMM_WORLD);
2684   t.start();
2685   my_function();
2686   MPI_Barrier (MPI_COMM_WORLD);
2687   t.stop();
2688   if (my MPI rank == 0)
2689     std::cout << "Calling my_function() took " << timer() << " seconds." << std::endl;
2690 ```
2691
2692 This way, you really only measure the time spent between when all
2693 processors have finished doing what they were doing before, and when they
2694 are all finished doing what they needed to do for `my_function`.
2695
2696 Another way to find some answers is to use the capabilities of the `Timer`
2697 class which can provide more detailed information when deal.II is
2698 configured to support MPI.
2699
2700 ## I have a special kind of equation!
2701
2702 ### Where do I start?
2703
2704 The deal.II tutorial has a number of programs that deal with particular
2705 kinds of equations, such as vector-valued problems, mixed discretizations,
2706 nonlinear and time-dependent problems, etc. The best way to start is to
2707 take a look at the existing tutorial programs and see if there is one that
2708 is already close to what you want to do. Then take that, try to understand
2709 its structure, and find a way to modify it to solve your problem as well.
2710 Most applications written based on deal.II are not written entirely from
2711 scratch, but have started out as modified tutorial programs.
2712
2713 ### Can I solve my particular problem?
2714
2715 The simple answer is: if it can be written as a PDE, then this is possible
2716 as evidenced by the many publications in widely disparate fields obtained
2717 with the help of deal.II. The more complicated answer is: deal.II is not a
2718 problem-solving environment, it is a toolbox that supports you in solving a
2719 PDE by the method of finite elements. You will have to implement assembling
2720 matrices and right hand side vectors yourself, as well as nonlinear outer
2721 iterations, etc. However, you will not need to care about programming a
2722 triangulation class that can handle locally refined grids in one, two, and
2723 three dimensions, linear algebra classes, linear solvers, different finite
2724 element classes, etc.
2725
2726 To give only a very brief overview of what is possible, here is a list of
2727 the nontrivial problems that were treated by the programs that the main
2728 authors alone wrote to date:
2729
2730  - Time-dependent acoustic and elastic wave equation, including nonlocal
2731    absorbing boundary conditions;
2732  - Stokes flow discretized with the discontinuous Galerkin finite element
2733    method;
2734  - General hyperbolic problems including Euler flow, using the
2735    discontinuous Galerkin finite element method;
2736  - Distributed parameter estimation problems;
2737  - Mixed finite element discretization of a mortar multiblock formulation
2738    of the Laplace equation;
2739  - Large-deformation elasto-plasticity in the simulation of plate
2740    tectonics.
2741
2742 To illustrate the complexity of the programs mentioned above we note that
2743 most of them include adaptive mesh refinement tailored to the efficient
2744 computation of specific quantities of physical interest and error
2745 estimation measured in terms of these quantities. This includes the
2746 solution of a so-called dual problems, that means e.g. for the wave
2747 equation the solution of a wave equation solved backward in time.
2748
2749 Problems other users of deal.II have solved include:
2750  - Elastoplasticity;
2751  - Porous media flow;
2752  - Crystal growth simulations;
2753  - Fuel cell simulations and optimization;
2754  - Fluid-structure interaction problems;
2755  - Time dependent large deformation problems for metal forming;
2756  - Contact problems;
2757  - Viscoelastic deformation of continental plates;
2758  - Glacial ice flows;
2759  - Thermoelastoplastic metal forming;
2760  - Eulerian coordinates problems in biomechanical modeling.
2761
2762 Some images from these applications can be found on this wiki's [[Gallery]]
2763 page. A good overview of the sort of problems that are being solved with the
2764 help of deal.II can also be obtained by looking at the large number of
2765 [publications](http://www.dealii.org/publications.html) written with the help of
2766 deal.II.
2767
2768 Probably, many other problem types are solved by the many users which we do
2769 not know of directly. If someone would like to have his project added to
2770 this page, just contact us.
2771
2772
2773 ### Why use deal.II instead of writing my application from scratch?
2774
2775 You can usually get the initial version of a code for any given problem
2776 done relatively quickly when you write it yourself, since the learning
2777 curve is not as steep as if you had to learn a new library; it's also true
2778 that it's easy to make this code twice as fast as if you had to use a
2779 library. In other words, this sounds like you should write finite element
2780 codes for your problem yourself.
2781
2782 However, you also need to keep in mind that it is the things you want to do
2783 after the first 3 months that will take you forever if you want to write it
2784 yourself, and where you will never be able to catch up with existing,
2785 established libraries: higher order elements; complicated, unstructured 3d
2786 meshes; parallelization; producing output in a format that's easily
2787 visualizable in 3d; adding an advected field for a tracer quantity; etc.
2788 Viewed this way, it's worth remembering that the primary commodity that's
2789 in short supply is not CPU time but your own programming time, and that's
2790 where you will be '''orders magnitude faster''' when using what others have
2791 already done, even if maybe your program ends up twice as slow as if you
2792 had written it from scratch with a particular application in mind.
2793
2794 ### Can I solve problems over complex numbers?
2795
2796 Yes, you can, and it has been done numerous times with deal.II. However, we
2797 have a standard recommendation: consider such problems as systems of
2798 partial differential equations, where the individual components of the
2799 solution are the real and imaginary part of your unknown. The reason for
2800 this is that for complex-valued problems, the product `<u,v>` of two vectors
2801 is not the same as `<v,u>`, and it is very easy to get this wrong in many
2802 places. If you want to avoid these common traps, then the easiest way
2803 around is to split up you equation into two equations of real and imaginary
2804 part first, and then treat the resulting system as a system of real
2805 variables. This also makes the type of linear system clearer that you get
2806 after discretization, and tells you something about which solver may be
2807 adequate for it.
2808
2809 The step-29 tutorial program shows how this is done for a complex-valued
2810 Helmholtz problem.
2811
2812 ### How can I solve a problem with a system of PDEs instead of a single equation?
2813
2814 The easiest way to do this is setting up a system finite element after you
2815 chose your base element, e.g.,
2816 ```cpp
2817 FE_Q<dim> base_element(2);
2818 FESystem<dim> system_element(base_element, 3);
2819 ```
2820
2821 will produce a biquadratic element for a system of 3 equations. With this
2822 finite element, all the functions that you always called for a scalar
2823 finite element should just work for this vector-valued element as well.
2824
2825 Refer to the step-8 and in particular to the step-20 tutorial programs for
2826 a lot more information on this topic. Several of the other tutorial
2827 programs beyond step-20 also use vector-valued elements and there is a
2828 whole module in the documentation on vector-valued problems that is worth
2829 reading.
2830
2831 ### Is it possible to use different models/equations on different parts of the domain?
2832
2833 Yes. The step-46 tutorial program shows how to do this: It solves a problem
2834 in which we solve Stokes flow in one part of the domain, and elasticity in
2835 the rest of the domain, and couple them on the interface. Similar
2836 techniques can be used if you want to exclude part of the domain from
2837 consideration, for example when considering voids in a body in which the
2838 governing equations do not make sense because there is no medium.
2839
2840 ### Can I solve problems with higher regularity requirements?
2841
2842 deal.II does not currently support *C<sup>1</sup>* or *C<sup>2</sup>* elements since these elements
2843 require specialized unit-to-reference mappings that have not yet been
2844 implemented. We recommend solving these problems by converting a higher-order
2845 problem into a lower-order one; e.g., one can rewrite a biharmonic problem as a
2846 coupled pair of Laplace equations.
2847
2848 Some projects have managed to use deal.II for problems with higher regularity
2849 requirements by exploiting properties of structured grids: see
2850 https://arxiv.org/abs/1810.02473 for more information.
2851
2852 ### Where do I start to implement a new Finite Element Class?
2853
2854 If you really need an element that isn't already implemented in deal.II,
2855 then you'll have to understand the interplay between FEValues, the finite
2856 element, the mapping, and quadrature objects. A good place to start would
2857 be to read the deal.II paper (Bangerth, Hartmann, Kanschat, ACM Trans.
2858 Math. Softw., 2007).
2859
2860 The actual implementation would most conveniently start from the `FE_Poly`
2861 class. You first implement the necessary polynomial space in the base
2862 library, then you derive `FE_Your_FE_Name` from `FE_Poly` (using your new
2863 polynomial class as a template) and add the connectivity information.
2864
2865 You'll probably need more specific help at various points -- this is what
2866 the mailing list is there for!
2867
2868 ## General finite element questions
2869
2870 ### How do I compute the error
2871
2872 If your goal is to compute the error in the form `||`u-u<sub>h</sub>`||` in some
2873 kind of norm, then you should use the function
2874 [VectorTools::integrate_difference](https://www.dealii.org/8.3.0/doxygen/deal.II/namespaceVectorTools.html#a01174a2a7e2ee8fa6abdfdd93ac7a317)
2875 which can compute the norm above in any number of norms (such as the L2, H1,
2876 etc., norms). Take a look at step-7.
2877
2878 On the other hand, if your goal is to *estimate* the error, then the one class
2879 that can do this is
2880 [Kelly Error Estimator](https://www.dealii.org/8.3.0/doxygen/deal.II/classKellyErrorEstimator.html).
2881 This class is used in most of the tutorial programs that use adaptively refined
2882 meshes, starting with step-6.
2883
2884 ### How to plot the error as a pointwise function
2885
2886 The functions mentioned in the previous question compute the error as a
2887 cellwise value. As a consequence, the values computed also include a factor
2888 that results from the size of the cell. If you're interested in the pointwise
2889 error as something that can be visualized, for example because you want to
2890 find a pattern in why the solution is not as you expect it to be, what you
2891 should do is this:
2892  - Interpolate the exact solution
2893  - Subtract the interpolated exact solution from the computed solution
2894  - Put the resulting vector into a
2895    [DataOut object](https://www.dealii.org/8.3.0/doxygen/deal.II/classDataOut.html).
2896    This will plot the nodal values of the errors u-u<sub>h</sub> on the current
2897    mesh.
2898
2899 As an example, the following code shows how to do this in principle:
2900 ```cpp
2901   template <int dim>
2902   class ExactSolution : public Function<dim>
2903   {
2904   public:
2905     ExactSolution () : Function<dim>(dim+1) {}
2906
2907     virtual double value (const Point<dim>   &p,
2908                           const unsigned int  component) const
2909     {
2910       return ...exact solution as a function of p...
2911     }
2912   };
2913
2914
2915   template <int dim>
2916   void MyProblem<dim>::plot_error () const
2917   {
2918     Vector<double> interpolated_exact_solution (dof_handler.n_dofs());
2919     VectorTools::interpolate (dof_handler,
2920                               ExactSolution<dim>(),
2921                               interpolated_exact_solution);
2922     interpolated_exact_solution -= solution;
2923
2924     DataOut<dim> data_out;
2925
2926     data_out.attach_dof_handler (dof_handler);
2927     data_out.add_data_vector (solution, "solution");
2928     data_out.add_data_vector (interpolated_exact_solution, "pointwise_error");
2929   }
2930 ```
2931
2932
2933 ### I'm trying to plot the right hand side vector but it doesn't seem to make sense!
2934
2935 In particular, what you probably see is that the plot shows values that are
2936 smaller by a factor of two along the boundary than in the inside, and by a
2937 factor of four in the corners (in 2d) or eight (in 3d). Similarly, on
2938 adaptively refined cells, the values appear to scale with the cell size.
2939 The reason is that trying to plot a right hand side vector doesn't make
2940 sense.
2941
2942 While you plot the vector as if it is function (by connecting dots with
2943 straight lines in 1d, or plotting surfaces in 2d), the thing you right hand
2944 side vector is in fact an element of the dual space. To wit:
2945  - A vector in primal space is a vector of nodal values so that `sum_i  U_i
2946    phi_i(x)` is a reasonable function of `x`. Solution vectors are examples
2947    of elements of primal space.
2948  - A vector in dual space is a vector W formed from the integration of an
2949    object in primal space against the shape functions, e.g. `W_i  = int
2950    f(x) phi_i(x)`. Examples of dual vectors are right hand side vectors.
2951
2952 For vectors in dual space, it doesn't make sense to plot them as functions
2953 of the form
2954      `sum_i W_i phi_i(x)`.
2955 The reason is that the values of the coefficients W_i are not of
2956 **amplitude** kind. Rather, the W_i are of kind amplitude (e.g. f(x)) times
2957 integration volume (the integral `*` dx over the support of shape functions
2958 phi_i). In other words, the sizes of cells comes into play for W_i, as does
2959 whether a shape function lies in the interior or at the boundary. In your
2960 case, the area of the integral when you integrate against shape functions
2961 at the boundary happens to be half the size of the integration area for
2962 shape functions in the interior.
2963
2964
2965 ### What does XXX mean?
2966
2967 The documentation of deal.II uses many finite element specific terms that
2968 may not always be entirely clear to someone not familiar with this
2969 language. In addition, we have certainly also invented our shares of
2970 deal.II specific terminology. If you encounter something you are not
2971 familiar with, take a look at the [deal.II glossary
2972 page](http://www.dealii.org/developer/doxygen/deal.II/DEALGlossary.html)
2973 that explains many of them.
2974
2975 ## I want to contribute to the development of deal.II!
2976
2977 deal.II is Open Source -- this not only implies that you as everyone else has
2978 access to the source codes, it also implies a certain development model:
2979 whoever would like to contribute to the further development is invited to do
2980 so: If you have changes or ideas, please send them to the
2981 [deal.II mailing list](http://www.dealii.org/mail.html)!
2982
2983 This model follows a small number of simple rules. The first and basic one
2984 is that if you have something that might be of interest to others as well,
2985 you are invited to send it to the list for possible inclusion into the
2986 library and use by others as well. Such additions useful to others are, for
2987 example:
2988  - new backends for output in a new graphical format;
2989  - input filters for some kind of data;
2990  - tool classes that do something that might be interesting to use in other
2991    programs as well.
2992
2993 A few projects (some easy, some difficult) can also be found in the
2994 [list of open issues](https://github.com/dealii/dealii/issues), where they
2995 are generally marked as "Enhancements".
2996
2997 If you consider providing some code for inclusion with the library, these
2998 are the simple rules of gaining reputation in the Open Source community:
2999  - your reputation grows with the number and complexity of your
3000    contributions;
3001  - your reputation with the maintainers of the library also grows with the
3002    degree of conformance of your proposed additions with the administrative
3003    rules stated below;
3004  - originators of code are credited full authorship.
3005
3006 In order to allow that a library remains a consistent piece of software,
3007 there are a number of administrative rules:
3008  - there are a number of maintainers that decide what goes into the
3009    library;
3010  - maintainers are benevolent, i.e. in general they want your addition to
3011    become part of the library;
3012  - however, they have to evaluate additions with respect to some criteria,
3013    among which are value for others;
3014  - whether it fits into the general framework (meaning that if your
3015    contribution requires the installation of some obscure other library
3016    that people do not usually have, then that must be discussed;
3017    alternatively, a way must be provided to disable your contribution on
3018    machines that do not have this lib);
3019  - completeness and amount of documentation;
3020  - existence and completeness of error checking through assertions.
3021
3022 However, again: the basic rule is that if you think your addition is
3023 interesting to others, there most probably is a way to get it into the
3024 library!
3025
3026
3027 ## I found a typo or a bug and fixed it on my machine. How do I get it included in deal.II?
3028
3029 First: thank you for wanting to do this! This software project is kept alive
3030 by people like you contributing to it. We like to include any improvement,
3031 even if it is just a single typo that you fixed.
3032
3033 If you have only a small change, or if this is your first time submitting
3034 changes, the easiest way to get them to us is by just emailing the
3035 [deal.II mailing lists](http://dealii.org/mail.html) and we will make sure they
3036 get incorporated. If you continue submitting patches (which we hope you will!)
3037 and become more experienced, we will start to ask you to use
3038 [git](http://en.wikipedia.org/wiki/Git_%28software%29) as the version control
3039 system and base your patches off of the
3040 [deal.II github repository](https://github.com/dealii/dealii).
3041
3042 The process for this is essentially the following (if you don't quite
3043 understand the terminology below, take a look at the manuals at the
3044 [github web site](https://github.com/), read
3045 [this online tutorial](https://www.atlassian.com/git/tutorial), or ask on the
3046 mailing list):
3047   - Create a github account
3048   - Fork the deal.II github repository, using the button at the top right
3049     of https://github.com/dealii/dealii
3050   - Clone the repository onto your local file system
3051   - Create a branch for your changes
3052   - Make your changes
3053   - Push your changes to your github repository
3054   - Create a pull request for your changes by going to your github
3055     account's deal.II tab where, after the previous step, there should be a
3056     button that allows you to create a pull request.
3057
3058 This list may sound intimidating at first, but in reality it's a fairly
3059 straightforward process that takes no more than 2 minutes after the first
3060 couple of times. But, as said, we'll be happy to hold your hand the first few
3061 times around and help you with the process! There's also a video lecture that
3062 demonstrates
3063 [how to submit a patch to github](http://www.math.colostate.edu/~bangerth/videos.676.32.8.html).
3064
3065 If you've submitted patches several times and know your way around git by now,
3066 please also consider to
3067   - make sure you base your patch off the most recent revision of the
3068     repository
3069   - you rewrite the history of your patch so that it contains a relatively
3070     small number of commits that are each internally consistent and could
3071     also be applied independently (see, for example, [the discussion
3072     towards the bottom of this
3073     page](https://github.com/dealii/dealii/pull/87)).
3074
3075
3076
3077 ## I'm fluent in deal.II, are there jobs for me?
3078
3079 Certainly. People with numerical skills are a sought commodity, both in
3080 academia and in businesses. In the US, the National Labs are also hiring
3081 lots of people in this field.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.