ENDIF()
#
-# As of 1.67, boost requires specifying the suffix for the python
-# component manually.
+# As of 1.67, boost requires specifying the suffix for the python
+# component manually.
#
UNSET(Boost_FOUND)
IF(${BOOST_VERSION} VERSION_LESS 1.67)
we add at least one new test to the testsuite, and we also do so if
we fix a bug, in both cases to make sure that future changes do not
break this functionality (again). Machines running the testsuite
- submit the results to our cdash
- <a href="https://cdash.43-1.org/index.php?project=deal.II">CDash instance</a>,
+ submit the results to our cdash
+ <a href="https://cdash.43-1.org/index.php?project=deal.II">CDash instance</a>,
a webpage showing the status of our regression tests.
</p>
# doxygen-generated HTML files.
SET(_last_year "1997")
FILE(GLOB _source_files
- "${CMAKE_SOURCE_DIR}/source/*/*.cc"
- "${CMAKE_SOURCE_DIR}/include/deal.II/*/*.h"
+ "${CMAKE_SOURCE_DIR}/source/*/*.cc"
+ "${CMAKE_SOURCE_DIR}/include/deal.II/*/*.h"
"${CMAKE_SOURCE_DIR}/examples/*/*.cc")
FOREACH (_source ${_source_files})
FILE(READ "${_source}" contents)
${_code_gallery_h_deps}
)
ADD_CUSTOM_TARGET(build_code-gallery_h
- DEPENDS
+ DEPENDS
${CMAKE_CURRENT_BINARY_DIR}/code-gallery.h
${_code_gallery_h_deps}
COMMENT
COMMAND ${PERL_EXECUTABLE}
ARGS
${CMAKE_SOURCE_DIR}/doc/doxygen/scripts/make_gallery.pl
- ${CMAKE_SOURCE_DIR}
- ${_step}
+ ${CMAKE_SOURCE_DIR}
+ ${_step}
${DEAL_II_CODE_GALLERY_DIRECTORY}/${_step}
- ${_relative_src_files}
+ ${_relative_src_files}
> ${CMAKE_CURRENT_BINARY_DIR}/${_step}.h
WORKING_DIRECTORY
${CMAKE_CURRENT_BINARY_DIR}
* A summary of the files that implement the interface to the supported auto-differentiable
* numbers is as follows:
*
- * - ad_drivers.h: Provides classes that act as drivers to the interface of internally supported
+ * - ad_drivers.h: Provides classes that act as drivers to the interface of internally supported
* automatic differentiation libraries. These are used internally as an intermediary to the
* helper classes that we provide.
* - ad_helpers.h: Provides a set of classes to help perform automatic differentiation in a
*
* @subsubsection auto_diff_1_3 User interface to the automatic differentiation libraries
*
- * The deal.II library offers a unified interface to the automatic differentiation libraries that
+ * The deal.II library offers a unified interface to the automatic differentiation libraries that
* we support. To date, the helper classes have been developed for the following contexts:
*
* - Classes designed to operate at the quadrature point level (or any general continuum point):
* - ResidualLinearization: Differentiation of a vector-valued finite element residual, leading to
* its consistent linearization.
*
- * Naturally, it is also possible for users to manage the initialization and derivative
+ * Naturally, it is also possible for users to manage the initialization and derivative
* computations themselves.
*
* The most up-to-date examples of how this is done using ADOL-C can be found in
* deal.II are now using features that were only introduced in C++11.
* That said, deal.II also has functions and classes that make using
* it with C++11 features easier.
- *
+ *
* One example is support for C++11
* <a href="http://en.wikipedia.org/wiki/C++11#Range-based_for_loop">range-based
* for loops</a>. deal.II-based codes often have many loops of the kind
* an introduction into the use of this class in the context of the
* vector-valued elasticity (Lamé) equations. step-20 discusses a mixed
* Laplace discretization that also uses vector-valued elements.
- *
+ *
* @ingroup feall
*/
* specify up front which information you are actually interested in. The
* UpdateFlags enumeration is used to offer symbolic names denoting what you
* want the FEValues class to compute.
- *
+ *
* All these classes are used in all
* @ref Tutorial "tutorial programs"
* from step-3 onward, and are described there in significant detail.
* <li> scalar: FE_DGP, FE_DGQ
* <li> scalar, different shape functions: FE_DGPMonomial, FE_DGPNonparametric, FE_DGQArbitraryNodes
* <li> vector-valued: FE_DGBDM, FE_DGNedelec, FE_DGRaviartThomas
- * </ul>
+ * </ul>
*
* @note The implementation of vector valued DG elements is supported
* by the class FE_DGVector, in the way, that only the vector
* it accepts a vector that describes a displacement field for each position
* of the domain. This is used in Eulerian computations without the need to
* actually move vertices after each time step.
- *
+ *
* In addition, the MappingC1 class provides for a boundary of the
* computational domain that is not only curved, but also has a continuous
* derivative at the interface between two cells on the boundary.
- *
+ *
* Finally, the MappingCartesian class is an optimization for elements that
* are brick-shaped and with edges parallel to the coordinate axes.
- *
+ *
* @ingroup feall
*/
/**
* @defgroup FE_vs_Mapping_vs_FEValues How Mapping, FiniteElement, and FEValues work together
- *
+ *
* <h2>Introduction</h2>
- *
+ *
* Most people create finite element (and, potentially, mapping) objects once
* but then never actually call any member functions on them -- they simply
* use them for assembly via the FEValues interface. The only other interaction
* values and derivatives on actual cells (obviously with the help of a
* mapping object). For the current discussion, only the latter role is
* important.
- *
+ *
* As with mappings, all that is important for us here is that the finite
* element classes can provide this information at given quadrature points,
* and that they can put the computed information into structures provided
* mapping class will indicate by adding update_contravariant_transformation
* to the list.
*
- *
+ *
* <h2>Pre-computing things</h2>
*
* At this point, the FEValues object has found out the complete
*
*
* <h3>Time dependent functions</h3>
- *
+ *
* For time dependent computations, boundary conditions and/or right
* hand side functions may also change with time. Since at a given
* time step one is usually only interested in the spatial dependence
* time of the function object to the present time step before handing
* it off to the VectorTools::interpolate_boundary_values method.
*
- *
+ *
* <h3>Tensor-valued functions</h3>
*
* The Function class is the most frequently used, but sometimes one needs a
correctness. Existing tutorial programs typically employ simpler rather than
more complicated solver schemes for exposition but frequently suggest more
complicated schemes including hints on how they might be implemented in an
- appendix.
+ appendix.
<li> <i>Try algorithms:</i> The rapid prototyping abilities of deal.II may
also help in determining best algorithms on the scale of programs to which
* allocates and manages degrees of freedom on vertices, faces, or cells, as
* described by the finite element class. This separation makes it possible to
* have multiple DoFHandler classes work on the same mesh at the same time.
- *
- *
+ *
+ *
* <h3>Grid generation</h3>
*
* There are three ways to create a mesh:
* The GridRefinement class implements a number of mesh refinement algorithms,
* based on refinement indicators given to its member functions.
*
- *
+ *
* <h3>Internal classes</h3>
*
* In addition to the above, there are a significant number of classes in this
* handling. They are generally in the internal namespace, and not meant for
* use in application code.
*
- *
+ *
* @author Wolfgang Bangerth, 1998-2006
*/
* for (unsigned int degree=1; degree<5; ++degree)
* fe_collection.push_back (FE_Q<dim>(degree));
* @endverbatim
- *
+ *
* This way, one can add elements of polynomial degree 1 through 4 to the
* collection. It is not necessary to retain the added object: the collection
* makes a copy of it, it does not only store a pointer to the given finite
* Instantiation of complex class and function templates is expensive both in
* terms of compile time and disk space. Therefore, we try to separate
* declaration and implementation of templates as far as possible, and make
- * sure that implementations are read by the compiler only when
+ * sure that implementations are read by the compiler only when
* necessary.
*
* Template classes in <tt>deal.II</tt> can be grouped into three categories,
* depending on the number of probable different instantiations. These three
* groups are discussed in the following.
- *
+ *
*
* @section Inst1 Known and fixed number of instantiations
*
* For these classes, adding instantiations for new parameters involves
* changing the library. However, this is rarely needed, of course, unless you
* are not content with computing only in 1d, 2d, or 3d.
- *
+ *
*
* @subsection Inst1a Available instances
*
* stored values: we only support <code>double</code>, <code>float</code>,
* and in some cases <code>std::complex@<double@></code> and
* <code>std::complex@<float@></code>.
- *
+ *
*
* @section Inst2 A few instantiations, most of which are known
*
* interfaces, including vector types that are not part of the library but
* possibly defined in an application program. In such a case, applications
* can instantiate these templates by hand as described in the next section.
- *
+ *
*
* @subsection Inst2c Creating new instances
*
* template class XXXX<long double>;
* @endcode
*
- *
+ *
* @subsection Inst2p Provided instances
*
* Like with the classes in section @ref Inst1, the instances provided in the
* itself), and it can create output used for visualization of data. Finally,
* it can output matrices in a graphical format.
*
- *
+ *
* <h3>Visualization of data</h3>
- *
+ *
* deal.II supports, through the DataOutBase class, a large number of popular
* visualization formats, such as those used by the OpenDX, gmv, or gnuplot
* programs. A complete list of supported formats is listed in the
* dependent simulations in the space-time domain: it collects the results
* from each time step and at the end outputs all of this information at once
* as a space-time file.
- *
- *
+ *
+ *
* <h3>Grid output</h3>
*
* Meshes, without any data vectors associated with it, can be written in a
*
* Finally, we have several auxiliary classes, namely MGLevelObject,
* which stores an object on each level*
- *
+ *
* See the step-16, step-16b, and step-39 example programs on how to use this
* functionality.
/**
* @defgroup PETScWrappers PETScWrappers
- *
+ *
* The classes in this module are
* wrappers around functionality provided by the PETSc library. They provide a
* modern object-oriented interface that is compatible with the interfaces of
/**
* Notations that reduce the order of tensors, effectively storing them in some
- * sort of consistent compressed storage pattern. An example is storing the
- * 6 independent components of $3\times 3$ symmetric tensors of rank 2 as a
- * vector with 6 components, and then representing the 36 independent elements
- * of symmetric $3\times 3 \times 3\times 3$ tensors of rank 4 (which when
- * applied to a symmetric rank-2 tensor yields another symmetric rank-2 tensor)
- * as a $6 \times 6$ matrix.
+ * sort of consistent compressed storage pattern. An example is storing the
+ * 6 independent components of $3\times 3$ symmetric tensors of rank 2 as a
+ * vector with 6 components, and then representing the 36 independent elements
+ * of symmetric $3\times 3 \times 3\times 3$ tensors of rank 4 (which when
+ * applied to a symmetric rank-2 tensor yields another symmetric rank-2 tensor)
+ * as a $6 \times 6$ matrix.
*
* Although this method of representing tensors is most regularly associated with
* the efficient storage of the fourth-order elasticity tensor, with its
- * generalization it has wider applicability. This representation is also common
+ * generalization it has wider applicability. This representation is also common
* in the physics, material science and FEM literature.
*
* There are several variations of tensor notation, each a slightly different
* structure. The primary difference between the various forms of tensor notation
* is the weighting prescribed to the various elements of the compressed tensors.
* This <a
- * href="https://en.wikipedia.org/wiki/Voigt_notation">wikipedia article</a> has
+ * href="https://en.wikipedia.org/wiki/Voigt_notation">wikipedia article</a> has
* some further general insights on this topic.
*
* @ingroup physics
* coordinate direction than in another) through the QAnisotropic class, as
* well as the definition of quadrature formulas that are not tensor products.
*
- *
+ *
* <h3>Use</h3>
- *
+ *
* Quadrature formulas are used, among other uses, when integrating matrix
* entries and the components of the right hand side vector. To this end, the
* quadrature point defined on the unit cell have to be mapped to the
* quadrature points as mapped to a real cell, for other uses as well. This
* can then be used, for example, to evaluate a right hand side function at
* these points.
- *
+ *
*
* <h3>QIterated</h3>
- *
+ *
* The class QIterated is used to construct an iterated quadrature formula out
* of an existing one, thereby increasing the accuracy of the formula without
* increasing the order. For example, by iterating the trapezoidal rule with
/**
* @defgroup SLEPcWrappers SLEPcWrappers
- *
+ *
* The classes in this module are wrappers around functionality
* provided by the SLEPc library. All classes and functions in this
* group reside in a namespace @p PETScWrappers.
/**
* @defgroup TrilinosWrappers TrilinosWrappers
- *
+ *
* The classes in this module are wrappers around functionality
* provided by the Trilinos library. They provide a modern
* object-oriented interface that is compatible with the interfaces of
/**
* @defgroup UpdateFlags The interplay of UpdateFlags, Mapping, and FiniteElement in FEValues
- *
+ *
* <h2>Introduction</h2>
- *
+ *
* In order to compute local contributions of an individual cell to the global
* matrix and right hand side, we usually employ two techniques:
* - First, the integral is transformed from the actual cell $K$ to the
* @f]
* where $q$ indicates the index of the quadrature point, $\hat{\bf x}_q$ its
* location on the reference cell, and $w_q$ its weight.
- *
+ *
* In order to evaluate such an expression in an application code, we
* have to access three different kinds of objects: a quadrature
* object that describes locations $\hat{\bf x}_q$ and weights $w_q$ of
* functions on the unit cell; and a mapping object that provides the
* Jacobian as well as its determinant. Dealing with all these
* objects would be cumbersome and error prone.
- *
+ *
* On the other hand, these three kinds of objects almost always appear together,
* and it is in fact very rare for deal.II application codes to do anything with
* quadrature, finite element, or mapping objects besides using them together.
* points and weights, mapped shape function values and derivatives as well as
* some properties of the transformation from the reference cell to the actual
* mesh cell.
- *
+ *
* Since computation of any of these values is potentially expensive (for
* example when using high order mappings with high order elements), the
* FEValues class only computes what it is explicitly asked for. To this
* property does not hold for the shape functions of Raviart-Thomas
* elements, however, which must be rotated with the local cell).
* This allows further optimization of the computations underlying assembly.
- *
+ *
*
* <h2> Tracking dependencies </h2>
- *
+ *
* Let's say you want to compute the Laplace matrix as shown above. In that
* case, you need to specify the <code>update_gradients</code> flag
* (for $\nabla\varphi_i(\bf x_q)$) and the <code>update_JxW_values</code>
* Jacobian matrix, $J^{-1}(\bf x_q)$ (and not just the determinant of the matrix),
* and to compute the inverse of the Jacobian, it is also necessary to compute
* the Jacobian matrix first.
- *
+ *
* Since these are requirements that are not important to the user, it
* is not necessary to specify this in user code. Rather, given a set
* of update flags, the FEValues object first asks the finite element
* store some temporary information that can be computed once and for
* all, and these flags will be used when re-computing data on each
* cell we will visit later on.
- *
- *
+ *
+ *
* <h2>Update once or each</h2>
- *
+ *
* As outlined above, we have now determined the final set of things that are
* necessary to satisfy a user's desired pieces of information as conveyed by
* the update flags they provided. This information will then typically be queried
* on every cell the user code visits in a subsequent integration loop.
- *
+ *
* Given that many of the things mappings or finite element classes compute are
* potentially expensive, FEValues employs a system whereby mappings and finite
* element objects are encouraged to pre-compute information that can be computed
* Jacobian of the mapping which depends on the geometry of the cell we visit;
* thus, for this element, simply copying pre-computed information is not
* sufficient to evaluate the values of shape functions on a particular cell.)
- *
+ *
* To accommodate this structure, both mappings and finite element classes
* may internally split the update flags into two sets commonly referenced as
* <code>update_once</code> and <code>update_each</code> (though these names
* set <code>update_once=update_values</code> and
* <code>update_each=0</code>, whereas the Raviart-Thomas element will
* do it the other way around.
- *
+ *
* These sets of flags are intended to be mutually exclusive. There is,
* on the other hand, nothing that ever provides this decomposition to
* anything outside the mapping or finite element classes -- it is a purely
* internal decomposition.
- *
- *
+ *
+ *
* <h2>Generation of the actual data</h2>
- *
+ *
* As outlined above, data is computed at two different times: once at
* the beginning on the reference cell, and once whenever we move to an
* actual cell. The functions involved in each of these steps are
* discussed next:
- *
- *
+ *
+ *
* <h3>Initialization</h3>
- *
+ *
* Computing data on the reference cell before we even visit the first
* real cell is a two-step process. First, the constructor of FEValues,
* FEFaceValues and FESubfaceValues, respectively, need to allow the
* <li>FiniteElement::get_face_data()
* <li>FiniteElement::get_subface_data()
* </ul>
- *
+ *
* The FEValues object then takes over ownership of these objects and will
* destroy them at the end of the FEValues object's lifetime. After this,
* the FEValues object asks the FiniteElement and Mapping objects to fill
* <li>FEFaceValues::initialize()
* <li>FESubfaceValues::initialize()
* </ul>
- *
- *
+ *
+ *
* <h3>Reinitialization for a mesh cell</h3>
- *
+ *
* Once initialization is over and we call FEValues::reinit, FEFaceValues::reinit
* or FESubfaceValues::reinit to move to a concrete cell or face, we need
* to calculate the "update_each" kinds of data. This is done in the following
* <li>FESubfaceValues::reinit() calls Mapping::fill_fe_subface_values(),
* then FiniteElement::fill_fe_subface_values()
* </ul>
- *
+ *
* This is where the actual data fields for FEValues, stored in
* internal::FEValues::MappingRelatedData and
* internal::FEValues::FiniteElementRelatedData objects, are
* computed. These functions call the function in Mapping first, such
* that all the mapping data required by the finite element is
* available. Then, the FiniteElement function is called.
- *
+ *
* @ingroup feall
*/
# "/a/b/c/name//doc//author", so make sure we eat the
# double slashes in the second step
SET_IF_EMPTY(DEAL_II_CODE_GALLERY_DIRECTORY ${CMAKE_SOURCE_DIR}/code-gallery)
-FILE(GLOB _code_gallery_names
+FILE(GLOB _code_gallery_names
"${DEAL_II_CODE_GALLERY_DIRECTORY}/*/doc/author")
STRING(REGEX REPLACE "/+doc/+author" "" _code_gallery_names "${_code_gallery_names}")
#
# Define target for the tutorial. It depends on the
# file tutorial.h built via the next target below, as well
-# as the various files we create from the tutorial
+# as the various files we create from the tutorial
# directories below that. These dependencies are added
# below the respective targets.
#
-<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
- "http://www.w3.org/TR/html4/loose.dtd">
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
+ "http://www.w3.org/TR/html4/loose.dtd">
-<html>
- <head>
+<html>
+ <head>
<title>The deal.II Readme on interfacing to OpenCASCADE</title>
<link href="../screen.css" rel="StyleSheet">
<meta name="copyright" content="Copyright (C) 2010 - 2015 by the deal.II authors">
- <meta name="keywords" content="deal.II">
- </head>
- <body>
+ <meta name="keywords" content="deal.II">
+ </head>
+ <body>
+
-
<h2>Installation of <acronym>OpenCASCADE</acronym></h2>
<p> <a href="http://www.opencascade.org/">OpenCASCADE</a> is a
<h5>BLAS and LAPACK</h5>
<p style="color: red">
- Note: At the time of writing (November 2016) Trilinos (at least up to
- v12.8.1) cannot link against LAPACK 3.6.0 or later, due to two symbols
- that were deprecated (and removed by default) in LAPACK 3.6.0 (see the
+ Note: At the time of writing (November 2016) Trilinos (at least up to
+ v12.8.1) cannot link against LAPACK 3.6.0 or later, due to two symbols
+ that were deprecated (and removed by default) in LAPACK 3.6.0 (see the
<a href="http://www.netlib.org/lapack/lapack-3.6.0.html">
- release notes</a>). To fix this, edit the Trilinos file
- <code>packages/epetra/src/Epetra_LAPACK_wrappers.h</code> and change the
+ release notes</a>). To fix this, edit the Trilinos file
+ <code>packages/epetra/src/Epetra_LAPACK_wrappers.h</code> and change the
lines
<pre>
-
+
#define DGGSVD_F77 F77_BLAS_MANGLE(dggsvd,DGGSVD)
#define SGGSVD_F77 F77_BLAS_MANGLE(sggsvd,SGGSVD)
</pre>
to
<pre>
-
+
#define DGGSVD_F77 F77_BLAS_MANGLE(dggsvd3,DGGSVD3)
#define SGGSVD_F77 F77_BLAS_MANGLE(sggsvd3,SGGSVD3)
</pre>
- before installing Trilinos. (Credit for this fix goes to
+ before installing Trilinos. (Credit for this fix goes to
<a href="https://github.com/gahansen/Albany/wiki/ALCF-Vesta">this page</a>.)
</p>
-
+
<p>
Trilinos sometimes searches for other libraries but can't find
them if they are not in the usual directories or have other
<a href="http://trilinos.org/docs/r11.8/packages/amesos/doc/html/index.html"
target="_top">this page for Trilinos 11.8</a>). Most of them are external
packages to Trilinos and you will need to tell Trilinos configuration
- scripts that you want to use them, for example via the
+ scripts that you want to use them, for example via the
<code>TrilinosWrappers::SolverDirect</code> class. This can be tricky,
but adding defines similar to the following to the cmake command line
will achieve the goal to enable the UMFPACK and SuperLU/SuperLUDist
refinement strategies based on local error indicators and error
estimators.</p>
-<li><p> Support for a variety of finite elements, including Lagrange
+<li><p> Support for a variety of finite elements, including Lagrange
elements of order one through four, and discontinuous elements.</p>
<li><p> Extensive documentation: all documentation is available online in a
- logical tree structure to allow fast access to the information you
+ logical tree structure to allow fast access to the information you
need. If printed it comprises about 200 pages of tutorials, several
reports, and far more than 1,000 pages of programming interface
- documentation with explanations of all classes, functions, and
+ documentation with explanations of all classes, functions, and
variables.</p>
<li><p> Modern software techniques that make access to the complex data
- structures and algorithms as transparent as possible. The use of
+ structures and algorithms as transparent as possible. The use of
object oriented programming allows for program structures similar to
the structures in mathematical analysis.</p>
<li><p> Fast algorithms that enable you to solve problems with up to several
- millions of degrees of freedom quickly. As opposed to programming
+ millions of degrees of freedom quickly. As opposed to programming
symbolic algebra packages the penalty for readability is low.</p>
-<li><p> Support for several output formats, including some common formats
+<li><p> Support for several output formats, including some common formats
for visualization of scientific data.</p>
<li><p> Support for a variety of computer platforms, including multi-
<li> <p>
Fix: getting the support points from system finite elements
(class <code>FESystem</code>) resulted in an
- exception.
+ exception.
<br>
(WB 2000/05/10)
</p>
</p>
<li> <p>
- Fix: slight bug in
+ Fix: slight bug in
<code>DataOut::build_patches</code>
in multithreaded mode fixed.
<br>
</p>
<li> <p>
- Fix:
+ Fix:
<code>SparsityPattern::print_gnuplot</code>
wrote rows and columns exchanged. Since most matrices have
symmetric sparsity patterns, this has gone unnoticed by now.
</p>
<li> <p>
- Fixed: the
- <code class="program">common/scripts/make_dependencies.pl</code>
+ Fixed: the
+ <code class="program">common/scripts/make_dependencies.pl</code>
script that sets up the dependencies for the make files had a
problem when the path to the library included special characters
such as `+'. This is now fixed.
</p>
<li> <p>
- New: In multithreading mode, the compile flags now also have
+ New: In multithreading mode, the compile flags now also have
<code class="program">_REENTRANT</code> defined. If this flag
is defined, then the standard Unix headers also declare
reentrant versions (with suffix <code class="program">_r</code>)
<li> <p>
New: the <code class="program">step-9</code> example program is
now ready, showing several advanced programming techniques in
- deal.II.
+ deal.II.
<br>
(WB 2000/07/18)
</p>
<li> <p>
- Changed: deal.II now uses the
+ Changed: deal.II now uses the
<code class="program">kdoc2</code> program to generate the API
documentation. This makes up for much nices documentation and
also works better than the previous <code class="program">kdoc1</code>.
</p>
<li> <p>
- Fixed: the
- <code class="program">common/scripts/make_dependencies.pl</code>
+ Fixed: the
+ <code class="program">common/scripts/make_dependencies.pl</code>
script that sets up the dependencies for the make files had a
problem when the path to the library included special characters
such as `+'. This is now fixed.
<li> <p>
Improved: the <code class="program">configure</code> script is
now used to select the compiler options. Previously, selection
- of compiler options was done both in
+ of compiler options was done both in
<code class="program">configure</code> as well as in the global
- options Makefile
+ options Makefile
<code class="program">common/Make.global_options</code>.
<br>
(WB 2000/06/02)
</p>
<li> <p>
- Improved: Dependence on the files
+ Improved: Dependence on the files
<code class="file">forward_declarations.h</code> has been widely
removed to improve compilation time
<br>
<li> <p>
- Fix: <code class="program">configure</code> now uses
+ Fix: <code class="program">configure</code> now uses
<code class="program">config.guess</code> to determine the
operating system and hardware platform.
<br>
<code>memory_consumption</code> that returns an
estimate of the amount of memory (in bytes) used up by this
class. Supporting functions to compute the size of STL vectors
- and other objects can be found the
+ and other objects can be found the
<code>MemoryConsumption</code> namespace.
<br>
(WB 2000/11/27)
faces in 3d computations).
<br>
(WB 2000/09/07)
- </p>
+ </p>
<li> <p>
New: There is a new <code>distance</code> function
distance between two points.
<br>
(Ralf Hartmann 2000/09/06)
- </p>
+ </p>
- <li> <p>
+ <li> <p>
New: <code>Timer</code> now uses the system
function <code>getrusage (RUSAGE_CHILDREN,
.)</code> that is need in multithreading. But still the <code
0 (at least on Solaris7).
<br>
(Ralf Hartmann 2000/08/25)
- </p>
+ </p>
<li> <p>
New: There are now a set of functions
product of tensors.
<br>
(WB 2000/07/23)
- </p>
+ </p>
<li> <p>
- New: The classes
+ New: The classes
<code>Patterns::Integer</code>
and <code>Patterns::Double</code> now allow that
a range may be specified in which the parameter shall
</p>
<li> <p>
- New: There are now functions
+ New: There are now functions
<code>SparsityPattern::symmetrize</code> and
<code>SparseMatrix::symmetrize</code> that
generate a symmetric matrix from a non-symmetric one.
(WB 2000/11/27)
</p>
- <li> <p>
+ <li> <p>
New: <code>
SparseMatrix<number>::n_actually_nonzero_elements</code>
returns the number of entries that are actually nonzero.
<br>
(Ralf Hartmann 2000/11/22)
- </p>
+ </p>
<li> <p>
Fixed: unlike announced in the docs, the
fixed now.
<br>
(GK 2000/08/28)
- </p>
+ </p>
<li> <p>
Improved: <code>PreconditionBlockSOR</code><code
of <code>SolverBicgstab::AdditionalData</code>.
<br>
(GK 2000/07/07)
- </p>
+ </p>
<li> <p>
New: <code>SolverRichardson</code> has a
(GK 2000/06/30)
</p>
- <li> <p>
+ <li> <p>
New: For recognizing a diverging solver before the maximum
- number of steps is reached,
- <code>SolverControl</code> returns
+ number of steps is reached,
+ <code>SolverControl</code> returns
<code>failure</code> also if the residual
increases over the start residual by a specific factor. This
- factor is given to the <code>SolverControl</code>
+ factor is given to the <code>SolverControl</code>
object by the <code>set_failure_criterion</code>
function. After calling the latter function, checking of this
additional failure criterion may again be disabled by calling
</p>
<li> <p>
- Fix:
+ Fix:
<code>SparsityPattern::print_gnuplot</code>
wrote rows and columns exchanged. Since most matrices have
symmetric sparsity patterns, this has gone unnoticed by now.
</p>
<li> <p>
- Fix: the
+ Fix: the
<code>FullMatrix::Tvmult</code> function
suffered from the same problems as the
<code>FullMatrix::Tmmult</code> function. This
<li> <p>
New: <code>BlockIndices</code>: Class that
manages the conversion of global indices into a block
- vector/matrix/... to the indices local to each of the blocks.
+ vector/matrix/... to the indices local to each of the blocks.
<br>
(WB 2000/05/08)
- </p>
+ </p>
<li> <p>
New: <code>BlockSparsityPattern</code> and
represent matrices that are composed of sparse matrices.
<br>
(WB 2000/05/08)
- </p>
+ </p>
<li> <p>
- Fix: the
+ Fix: the
<code>FullMatrix::mmult</code> and
<code>FullMatrix::Tmmult</code> code don't
resize their output argument any more, as this is not common
- style in the library. Furthermore,
+ style in the library. Furthermore,
<code>FullMatrix::Tmmult</code> was utterly
- broken.
+ broken.
<br>
(WB 2000/05/08)
</p>
<li> <p>
Change: the <code>matrix_norm</code> functions
- of sparse and full matrices are renamed to
+ of sparse and full matrices are renamed to
<code>matrix_norm_square</code>, since they in
fact return the square of the norm. This should avoid confusion
in some cases.
</p>
<li> <p>
- Fix: the ``copy-like'' constructor of
+ Fix: the ``copy-like'' constructor of
<code>SparsityPattern</code> that copies another
object and adds some off-diagonals had a bug that caused an
exception in some cases. This is now fixed.
</p>
<li> <p>
- New: there are now functions <code>Vector::swap</code>
+ New: there are now functions <code>Vector::swap</code>
and <code>BlockVector::swap</code>, as well as
global functions <code>swap(u,v)</code> that
exchange the data of two vectors without needing a temporary
<li> <p>
Change: the solver classes in LAC lost their first template
- argument. Their names are now
+ argument. Their names are now
<code>SolverXX<VECTOR></code>, where
<code>XX</code> denotes the name of the solver
(e.g. CG, GMRES, etc). Furthermore, the
<li> <p>
Extend: <code>DoFTools::extract_boundary_dofs</code>
now allows to also specify which boundary conditions shall be
- considered.
+ considered.
<br>
(WB 2000/12/04)
</p>
<li> <p>
- New: some arguments of the functions
+ New: some arguments of the functions
<code>DoFHandler::n_boundary_dofs</code>,
<code>DoFTools::extract_boundary_dofs</code>,
and
- <code>DoFTools::map_dof_to_boundary_indices</code>
+ <code>DoFTools::map_dof_to_boundary_indices</code>
are changed from <code>list</code> to
<code>set</code>, since that resembles more
closely the purpose of the parameter, and makes computations
class can handle boundary information in 3d as well
(i.e. lines and quads in 3d with special material IDs).
<br>
- (<a href="mailto://ms@biomech.tu-graz.ac.at">Michael
+ (<a href="mailto://ms@biomech.tu-graz.ac.at">Michael
Stadler</a> 2000/10/11)
- </p>
+ </p>
<li> <p>
Extended: The <code>GridIn</code>
<br>
(<a href="mailto://ms@biomech.tu-graz.ac.at">Michael
Stadler</a> 2000/10/11)
- </p>
+ </p>
<li> <p>
New: The <code>GridRefinement::refine_and_coarsen_optimize</code>
refinement and coarsening.
<br>
(Thomas Richter 2000/10/10)
- </p>
+ </p>
<li> <p>
Extended: The <code>GridIn::delete_unused_vertices</code>
to the triangulation object.
<br>
(WB 2000/09/26)
- </p>
+ </p>
<li> <p>
New: The <code>GridIn</code>
class can now read the basics of grids in DB Mesh format.
<br>
(WB 2000/09/26)
- </p>
+ </p>
<li> <p>
Extended: The <code>KellyErrorEstimator</code>
<code>DoFHandler</code> object.
<br>
(WB 2000/09/11)
- </p>
+ </p>
<li> <p>
New: The <code>DataOut_Faces</code> class allows to
cuts through the domain.
<br>
(WB 2000/09/07)
- </p>
+ </p>
<li> <p>
- Removed: The
+ Removed: The
<code>DataOut_Old</code> class has finally gone for
good. It was already deprecated in version 3.0, and has been
superceded for a long time by the framework of classes around
</p>
<li> <p>
- Changed: The
+ Changed: The
<code>map_dof_to_boundary_index</code>
functions have been moved from the <code>DoFHandler</code>
to the <code>DoFTools</code> class, in order to
were dormant, since the values of finite elements are always
computed for the presently available Lagrange elements, but
would have been activated once there are other classes of
- elements.
+ elements.
<br>
(Ralf Hartmann 2000/07/20)
</p>
<br>
(GK 2000/07/12)
</p>
-
+
<li> <p>
New: The new function <code>get_face</code>
of <code>FEFaceValues</code> and <code
</p>
<li> <p>
- Fix: slight bug in
+ Fix: slight bug in
<code>DataOut::build_patches</code>
in multithreaded mode fixed.
<br>
</p>
<p>
- Removed: The obsolete
+ Removed: The obsolete
<code>MatrixCreator::create_interpolation_matrix</code>
function is now removed.
</p>
</p>
<li> <p>
- Changed: enum
- <code>MeshSmoothing</code> is moved into the
+ Changed: enum
+ <code>MeshSmoothing</code> is moved into the
<code>Triangulation</code> class.
<br>
(Ralf Hartmann 2000/05/18)
</p>
- <li> <p>
+ <li> <p>
New: <code>Triangulation<dim>::patch_level_1</code>
is a new mesh smoothing. A mesh of patch level 1 consists of
patches, i.e. they consists of cells that are all refined at
</p>
<li> <p>
- Changed: the
+ Changed: the
<code>MatrixTools::apply_boundary_values</code>
now uses a much faster algorithm when working on matrices with
symmetric sparsity patterns. On the other hand, it does no more
eliminate whole rows when a matrix has a non-symmetric sparsity
pattern, or if the user (through a new flag) tells the function
that this is not necessary, for example if the matrix itself is
- non-symmetric.
+ non-symmetric.
</p>
<p>
</p>
<li> <p>
- Changed: The computation of the Jacobian matrices in the
+ Changed: The computation of the Jacobian matrices in the
<code>FEValues</code> class is now done more
efficiently. The speedup is in the range of a factor of 40 for
- 3D.
+ 3D.
<br>
(John Burnell, WB 2000/05/16)
</p>
- <li> <p>
+ <li> <p>
Change: <code>DoFTools::make_hanging_node_constraints ()</code>
does not use the user flags any more, and can thus run in
- parallel more than once.
+ parallel more than once.
<br>
(WB 2000/05/15)
</p>
- <li> <p>
+ <li> <p>
Extended: <code>DoFTools::make_sparsity_pattern ()</code>
now accepts a template parameter as sparsity pattern. This
- allows to use this function for the usual
- <code>SparsityPattern</code>, or for
+ allows to use this function for the usual
+ <code>SparsityPattern</code>, or for
<code>BlockSparsityPattern</code> arguments.
<br>
(WB 2000/05/15)
</p>
- <li> <p>
+ <li> <p>
New: <code>DoFTools::extract_hanging_node_dofs ()</code>
identifies nodes that will be constrained by hanging node constraints.
<br>
(WB 2000/05/12)
</p>
- <li> <p>
+ <li> <p>
New: <code>DoFRenumbering::sort_selected_dofs_back ()</code>
sorts selected degrees of freedom to the end of the index
range.
(WB 2000/05/12)
</p>
- <li> <p>
+ <li> <p>
Change: the return value of
<code>DoFHandler::max_couplings_between_dofs ()</code>
is bounded by <code>DoFHandler::n_dofs()</code>.
(Ralf Hartmann 2000/05/11)
</p>
- <li> <p>
+ <li> <p>
New: <code>FEValuesBase::get_cell ()</code>
returns present cell.
<br>
(Ralf Hartmann 2000/05/11)
</p>
- <li> <p>
+ <li> <p>
Fix: <code>FESystem::reinit()</code> generated
- an exception if <code>update_support_points</code>
+ an exception if <code>update_support_points</code>
was set.
<br>
(WB 2000/05/10)
</p>
- <li> <p>
- New: <code>IntergridMap::get_{source,destination}_grid</code>
+ <li> <p>
+ New: <code>IntergridMap::get_{source,destination}_grid</code>
functions return the grids for which the map was created.
<br>
(WB 2000/04/19)
</p>
<li> <p>
- New: <code>DoFTools::extract_boundary_dofs</code>
+ New: <code>DoFTools::extract_boundary_dofs</code>
finds all degrees of freedom which are at the boundary and belong to
specified components.
<br>
</p>
<li> <p>
- New: <code>DoFTools::compute_intergrid_constraints</code>
+ New: <code>DoFTools::compute_intergrid_constraints</code>
allows to use different discretization grids for different
variables.
<br>
<li> <p>
New: <code>DataOut::clear_data_vectors</code>
allows to re-use an object without deleting the <code
- class="class">DoFHandler</code>.
+ class="class">DoFHandler</code>.
<br>
(GK 2000/04/05)
</p>
<br>
(<a href="mailto:or@winfos.com">Oliver Rheinbach</a> 2001/04/12)
</p>
-
+
<li> <p>
- Fixed: When using Neuman boundary functions in the
+ Fixed: When using Neuman boundary functions in the
<code>KellyErrorEstimator</code> class, it was
assumed that the function object had <code
class="member">Function::vector_value</code> overloaded, even
<br>
(WB 2001/04/09)
</p>
-
+
<li> <p>
New/Fixed: Now there exists a new <code
class="class">Triangulation</code>::<code
<br>
(RH 2001/04/03)
</p>
-
+
<li> <p>
Fixed: Missing <code>ucd</code> is now added to the list of
supported output formats returned by <code
(Roger Young, WB 2001/03/22)
</p>
- <li> <p>
+ <li> <p>
Removed: The explicit instantiations of <code
class="class">SparseMatrix<long double></code> are removed as a
prerelease of gcc3.0 fails to compile it. Now the user of <code
Changed: We now unconditionally include
<code>deal.II/include/grid/tria_accessor.templates.h</code>
(which contains some inline functions for triangulation
- accessor classes) into
+ accessor classes) into
<code>deal.II/include/grid/tria_accessor.h</code> to work
around a problem with gcc3.0 which does not place a copy of
these functions into the library. Previously we only included
<li> <a name="new_fe_mapping_design"></a>
<i> New Design of <code>FiniteElements</code>
and <code>Mappings</code></i>
-
+
<p>
New: The finite element
classes have been redesigned from scratch to allow also
to higher degree.
</p>
- <p>
+ <p>
Furthermore we have totally reimplemented the mapping
between unit and real cell. In future we won't be restricted
to Q<sub>1</sub> mapping any more, but we will have
</p>
<li> <p>
- New: There is now some support to include and use routines from the
- <a href="http://www.cse.clrc.ac.uk/Activity/HSL"
+ New: There is now some support to include and use routines from the
+ <a href="http://www.cse.clrc.ac.uk/Activity/HSL"
target="_top">Harwell Subroutine Library</a>.
<br>
(WB 2001/01/30)
<h3>base</h3>
<ol>
- <li> <p>
+ <li> <p>
New: Color function <code
class="member">DataOutBase::EpsFlags::reverse_grey_scale_color_function</code>.
<br>
(WB 2001/08/24)
</p>
- <li> <p>
+ <li> <p>
New: Function <code>QProjector::project_to_child</code>
generates quadrature formulae which act on the area which a
child would occupy.
(WB 2001/07/18)
</p>
- <li> <p>
+ <li> <p>
New: classes <code>FourierSineFunction</code>
and <code>FourierCosineFunction</code>,
resembling one mode of a Fourier decomposition. Classes <code
New: Classes <code>Vector</code> and
<code>BlockVector</code> now have member functions
<code>operator*=</code> which scale the vectors
- by a constant factor. These functions are aliases for the
+ by a constant factor. These functions are aliases for the
<code>scale</code> member functions except that
they return a reference to themselves.
<br>
</p>
<li> <p>
- New: There is now a function
- <code>FullMatrix::precondition_Jacobi</code>. The
+ New: There is now a function
+ <code>FullMatrix::precondition_Jacobi</code>. The
<code>PreconditionJacobi</code> class is
therefore now also applicable with the full matrix class.
<br>
</p>
<li> <p>
- New: There are now functions
+ New: There are now functions
<code>SparsityPattern::copy_from</code> and
<code>SparseMatrix::copy_from</code>
that can be used to construct sparsity patterns and matrix at
once, i.e. without setting up a sparsity pattern element by
element, possibly after using a way too large value for the
- maximal row length, then calling
+ maximal row length, then calling
<code>SparsityPattern::compress</code> etc.
<br>
(WB 2001/05/07)
(GK 2001/03/29)
</p>
- <li> <p>
+ <li> <p>
New: The functions <code>FullMatrix</code>::<code
class="member">mmult</code> and <code
class="member">Tmmult</code> now have an additional
New: class <code>SchurMatrix</code> implements
a Schur complement for block matrices. It provides matrix-vector
multiplication suitable for iterative methods as well as pre- and
- post-processing of right hand side and slution, respectively.
+ post-processing of right hand side and slution, respectively.
<br>
(GK 2001/03/22)
</p>
<br>
(RH 2001/03/14)
</p>
-
- <li> <p>
+
+ <li> <p>
New: class <code
class="class">BlockMatrixArray<MATRIX></code> implements
a block matrix based on an array of matrix pointers. Since this
</p>
<li> <p>
- New: There is now some support to include and use routines from the
- <a href="http://www.cse.clrc.ac.uk/Activity/HSL"
+ New: There is now some support to include and use routines from the
+ <a href="http://www.cse.clrc.ac.uk/Activity/HSL"
target="_top">Harwell Subroutine Library</a>, and support
- classes
+ classes
<code>SparseDirectMA27</code> and
<code>SparseDirectMA47</code>
for the sparse direct solvers MA27 and MA47.
<h3>deal.II</h3>
<ol>
- <li> <p>
+ <li> <p>
New: Class <code>MappingQ1Eulerian</code>
implementing an Eulerian mapping.
<br>
(<a href="mailto:ms@biomech.tu-graz.ac.at">Michael Stadler</a> 2001/09/24)
</p>
- <li> <p>
+ <li> <p>
New: <code>VectorTools::create_boundary_right_hand_side</code>
integrates boundary forces for inhomogeneous Neumann boundary values.
<br>
(WB 2001/09/13)
</p>
- <li> <p>
+ <li> <p>
New: <code>DoFTools::make_flux_sparsity_pattern</code>
now exists also for 1d.
<br>
(WB 2001/09/03)
</p>
- <li> <p>
+ <li> <p>
New: There are now two functions
<code
class="member">FETools::hierarchic_to_lexicographic_numbering</code>
(WB 2001/08/31)
</p>
- <li> <p>
+ <li> <p>
New: <code>ConstraintMatrix::close</code>
now simply returns instead of throwing an exception, if the
matrix was already closed.
New: Implement
<code>DoFTools::make_sparsity_pattern</code>,
<code>DoFTools::make_boundary_sparsity_pattern</code>,
- and
+ and
<code>ConstraintMatrix::condense</code> to work on
the <code>CompressedSparsityPattern</code> and
<code>CompressedBlockSparsityPattern</code> and
<li> <p>
New: There are now functions <code
class="class">DoFTools</code>::<code
- class="member">map_dofs_to_support_points</code> and
+ class="member">map_dofs_to_support_points</code> and
<code
class="class">DoFTools</code>::<code
class="member">map_support_points_to_dofs</code> that
<br>
(<a href="mailto:roger@kea.grace.cri.nz">Roger Young</a> 2001/04/12)
</p>
-
+
<li> <p>
Fixed: the <code
class="class">DoFRenumbering</code>::<code
<br>
(<a href="mailto:or@winfos.com">Oliver Rheinbach</a> 2001/04/12)
</p>
-
+
<li> <p>
Changed: A <code>typedef FunctionMap</code> was
declared in at least four places previously (in classes <code
<br>
(WB 2001/04/09)
</p>
-
+
<li> <p>
- Fixed: When using Neuman boundary functions in the
+ Fixed: When using Neuman boundary functions in the
<code>KellyErrorEstimator</code> class, it was
assumed that the function object had <code
class="class">Function</code>::<code
<br>
(WB 2001/04/09)
</p>
-
+
<li> <p>
New: The various <code
class="class">MatrixCreator</code>::<code
<br>
(WB 2001/04/08)
</p>
-
+
<li> <p>
New: The following functions are extended to work
on arbitrary mappings:
<p>
All these functions now take a <code
class="class">Mapping</code> object as additional argument.
- <br>
+ <br>
For backward compatibility there exist second versions of
each of these functions that still don't have a <code
class="class">Mapping</code> argument. These functions
<br>
(RH 2001/04/04)
</p>
-
+
<li> <p>
New/Fixed: Now there exists a new <code
class="class">Triangulation</code>::<code
<br>
(RH 2001/04/03)
</p>
-
+
<li> <p>
Fixed: Missing <code>ucd</code> is now added to the list of
supported output formats returned by <code
<li> <p>
Removed: the <code>ProblemBase</code> class,
- which has been deprecated since before the release of
+ which has been deprecated since before the release of
deal.II 3.0, has finally been removed. The
same applies for the classes
<code>Assembler</code>,
the class <code>VectorTools</code> are extended
to work on arbitrary mappings. All these functions now take a
<code>Mapping</code> object as additional
- argument.
+ argument.
<br>
For backward compatibility there exist second versions of
each of these functions that still don't have a <code
(WB 2001/03/14)
</p>
- <li> <p>
+ <li> <p>
Removed: The function <code
class="class">FiniteElement<dim></code>::<code
class="member">get_support_points</code>
(RH 2001/03/14)
</p>
- <li> <p>
+ <li> <p>
New: The class <code>Boundary</code> has two new
functions <code
class="member">get_intermediate_points_on_line</code> and <code
<br>
(RH 2001/03/14)
</p>
-
+
<li> <p>
Changed: We now unconditionally include
<code>deal.II/include/grid/tria_accessor.templates.h</code>
(which contains some inline functions for triangulation
- accessor classes) into
+ accessor classes) into
<code>deal.II/include/grid/tria_accessor.h</code> to work
around a problem with gcc3.0 which does not place a copy of
these functions into the library. Previously we only included
(WB 2001/02/26)
</p>
- <li> <p>
+ <li> <p>
New: The <code>DoFTools</code> class now has
a function <code>count_dofs_per_component</code>
that counts the number of degrees of freedom in each of the
<ol>
<li> <p>
- New: Output for
+ New: Output for
<a href="http://www.amtec.org" target="_top">Tecplot</a> has
been added. It can be used by choosing output format «tecplot».
<br>
</p>
<li> <p>
- New: Output for
+ New: Output for
<a href="http://www.opendx.org" target="_top">OpenDX</a> has
been added. It can be used by choosing output format «dx» (not
yet for grid output). The data format is very basic now, but it
<li> <p>
New: the <code>SparseDirectMA27/47</code> classes
now provide access to Mutex locks for external
- synchronisation.
+ synchronisation.
<br>
(WB 2001/11/14)
</p>
class="member">GridRefinement::refine_fixed_fraction</code>
function sometimes had problems when indicators vary over
several orders of magnitude, due to roundoff. This is now
- fixed.
+ fixed.
<br>
(WB 2001/11/05)
</p>
- <li> <p>
+ <li> <p>
New: <code
class="member">DoFTools::extract_subdomain_dofs</code> selects
those degrees of freedom which are located on cells with a
<ol>
- <li> <p>
+ <li> <p>
Fixed: In 3d, the function <code
- class="member">DoFTools::make_hanging_node_constraints</code>
+ class="member">DoFTools::make_hanging_node_constraints</code>
contained an assertion that failed erroneously for finite
elements that do not have degrees of freedom on vertices. This
is now fixed.
- <br>
+ <br>
(WB 2002/02/21)
</p>
<li> <p>
Fixed: <code>TriaAccessor<3,3>::measure</code>
sometimes computed a negative value. This is now fixed.
- <br>
+ <br>
(WB 2002/02/21)
</p>
</ol>
(BK 2002/05/22)
</p>
- <li> <p>
+ <li> <p>
New: The step-14 example program demonstrates duality based
error estimators and some more software design concepts.
<br>
(WB 2002/05/05)
</p>
- <li> <p>
+ <li> <p>
New: In all previous versions, deal.II used
the <a href="http://www.cs.wustl.edu/~schmidt/ACE.html"
target="_top">ACE (Adaptive Communications Environment)</a>
available on many systems. The relieves you from the need of
installing a huge library of which the most part is not used
anyway. However, if you use ACE for other reasons, then it is
- still supported. For installation instructions, see the
+ still supported. For installation instructions, see the
<a href="../../readme.html" target="body">ReadMe</a> file.
<br>
(WB 2002/04/30)
</p>
- <li> <p>
+ <li> <p>
Changed: The Makefiles for the library are now truly
parallel. To this end, the automatic generation of the files
<tt>forward_declarations.h</tt> in the various directories had
(GK 2002/04/17)
</p>
- <li> <p>
+ <li> <p>
New: The step-13 example program tells you something about
software design things for finite element programs.
<br>
(WB 2002/04/16)
</p>
- <li> <p>
+ <li> <p>
Changed: Due to problems with undeclared functions and general
compatibility concerns, <code>-ansi</code> is now no more part
of the compile flags for gcc.
(WB 2002/04/16)
</p>
- <li> <p>
+ <li> <p>
Fixed: Explicit specializations of member templates are now
conforming to the C++ standard. While most compilers accepted
the previous form, Sun's Forte compiler wants a strictly
(WB 2002/03/25)
</p>
- <li> <p>
+ <li> <p>
Fixed: For gcc versions that used <code
class="class">ostrstream</code> instead of <code
class="class">ostringstream</code>, it was necessary to append
(WB 2002/03/13)
</p>
- <li> <p>
+ <li> <p>
Changed: The configure machinery has been revamped
significantly.
<br>
(WB 2002/03/08)
</p>
- <li> <p>
+ <li> <p>
Added: The top-level Makefile now supports "optimized" as a
target that builds only optimized versions of the <code>base</code>,
<code>lac</code>, <code>1d</code>, <code>2d</code>, and <code>3d</code>
- libraries.
+ libraries.
<br>
(BK 2002/02/19)
</p>
- <li> <p>
+ <li> <p>
Changed: The build system was entirely revised. Object
files in debug mode now have the suffix <code>.g.o</code>
instead of <code>.go</code>. All object files from the
<h3>base</h3>
<ol>
- <li> <p>
+ <li> <p>
New: The <code>vector2d</code> row accessor
classes now have member functions <code
class="member">begin</code> and <code>end</code>
which allow iterating over the elements of a row of such an
- object.
+ object.
<br>
(WB 2002/05/30)
</p>
- <li> <p>
+ <li> <p>
New: The <code>Legendre</code> and
<code>LagrangeEquidistant</code> classes now have
static member functions <code
(WB 2002/05/27)
</p>
- <li> <p>
+ <li> <p>
Changed: The <code>Polynomial</code> and
<code>LagrangeEquidistant</code> classes have lost
their default constructor, as that did not make much sense
(WB 2002/05/27)
</p>
- <li> <p>
+ <li> <p>
Fixed: When forward declaring the <code
class="class">Tensor</code> class, we now also forward declare
its partial specialization for a rank one tensor. Not doing so
(WB 2002/03/22)
</p>
- <li> <p>
+ <li> <p>
Fixed: The class <code>TensorFunction</code>
now uses local types <code>value_type</code> and
<code>gradient_type</code> as return values of
(WB 2002/03/20)
</p>
- <li> <p>
+ <li> <p>
Improved: The <code>AssertThrow</code> macro now
uses <code>__builtin_expect</code> if the
compiler supports this. This indicates to the compiler that we
(i.e. like to a two-dimensional C-style array).
<br>
(WB 2002/03/08)
- </p>
+ </p>
<li> <p>
Changed: The function <code>MultithreadInfo</code>::
of CPUs when running on Silicon Graphics.
<br>
(BK 2002/02/19)
- </p>
+ </p>
- <li> <p>
+ <li> <p>
Changed: The quite logorrhoeic function name <code
class="class">TensorProductPolynomials</code>::<code
class="member">n_tensor_product_polynomials</code> was changed to
(GK 2002/02/11)
</p>
- <li> <p>
+ <li> <p>
New: The class <code>PolynomialSpace</code>
implements the space of polynomials at most a certain degree in
arbitrary space dimensions.
(GK 2002/02/11)
</p>
- <li> <p>
+ <li> <p>
New: The function <code>DataOutBase</code>::
<code>write_tecplot_binary</code> has been
added. This function will write Tecplot binary files if the
(GK 2002/05/24)
</p>
- <li> <p>
+ <li> <p>
New: Function <code
class="member">SparsityPattern::matrix_position</code> is the
inverse function for <code
(WB 2002/02/13)
</p>
- <li> <p>
+ <li> <p>
New: Functions <code
class="member">SparsityPattern::copy_from</code> and <code
class="member">SparseMatrix::copy_from</code> allow to copy a full
<h3>deal.II</h3>
<ol>
- <li> <p>
+ <li> <p>
New: The <code>GeometryInfo</code> class now
provides two methods,
<code>unit_cell_vertex</code> and <code
(GK 2002/05/29)
</p>
- <li> <p>
+ <li> <p>
New: The <code>GridOut::</code>
<code>write_dx</code> function is now implemented.
It allows to write the mesh (cells and faces) with some additional
(GK 2002/05/02)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code>IteratorState::IteratorState</code>
enum is now called <code
class="class">IteratorState::IteratorStates</code>. This works
(WB 2002/03/20)
</p>
- <li> <p>
+ <li> <p>
Fixed: In 3d, the function <code
- class="member">DoFTools::make_hanging_node_constraints</code>
+ class="member">DoFTools::make_hanging_node_constraints</code>
contained an assertion that failed erroneously for finite
elements that do not have degrees of freedom on vertices. This
is now fixed.
- <br>
+ <br>
(WB 2002/02/21)
</p>
- <li> <p>
+ <li> <p>
Fixed: <code>TriaAccessor<3,3>::measure</code>
sometimes computed a negative value. This is now fixed.
- <br>
+ <br>
(WB 2002/02/21)
</p>
- <li> <p>
+ <li> <p>
New: Finite element family with complete polynomial spaces
for discontinuous Galerkin: <code>FE_DGP</code>
- <br>
+ <br>
(GK 2002/02/11)
</p>
</ol>
<h3>General</h3>
<ol>
- <li> <p>
+ <li> <p>
New: deal.II now uses a new threading
scheme. The new scheme is simpler to use, and in particular
more flexible in some cases where only one thread is started,
(WB 2003/02/06)
</p>
- <li> <p>
+ <li> <p>
New: deal.II now makes use of some parts of
the <a href="http://www.boost.org/">boost</a> library, which is
supposed to be a testground for the next generation C++ standard
the string <code>_local</code> since you may still want to use
another version or installation of boost in your own programs.
</p>
-
+
<p>
Also note that boost is large -- much larger than the subset we
have imported --, so we only took what we needed.
(WB 2003/02/06)
</p>
- <li> <p>
+ <li> <p>
Fixed: A longstanding bug in the documentation system has been fixed: if
a namespace was ended with <code>}</code> instead of <code>};</code>
(note the semicolon), then the documentation tool assumed that the
(WB 2003/02/02)
</p>
- <li> <p>
+ <li> <p>
Changed: Classes and structures that are declared inside namespaces named
<code>internal</code> are now no longer shown in the class
index view of the documentation of each of the sublibraries. Since they
(WB 2003/02/02)
</p>
- <li> <p>
- Fixed: Some of the formulas in the step-14 tutorial were obviously
+ <li> <p>
+ Fixed: Some of the formulas in the step-14 tutorial were obviously
scrambled a little. This is now fixed.
<br>
(Roy Stogner, WB 2003/01/30)
</p>
- <li> <p>
- Changed: The main Makefile has been changed to sequentialize building
+ <li> <p>
+ Changed: The main Makefile has been changed to sequentialize building
the base, lac, and deal.II sublibraries. We changed this, since on some
systems (notably AIX), the latter libraries need to be linked against
the former ones, when creating shared libraries.
(WB 2003/01/24)
</p>
- <li> <p>
- New: Changes have been made to support compiling and using
+ <li> <p>
+ New: Changes have been made to support compiling and using
deal.II on AIX 5 systems.
<br>
(WB 2003/01/24)
</p>
- <li> <p>
+ <li> <p>
Removed: Thread support now relies solely on the use of POSIX
functions. The use of the
<a href="http://www.cs.wustl.edu/~schmidt/ACE.html" target="_top">ACE
(WB 2003/01/11)
</p>
- <li> <p>
+ <li> <p>
New: Some changes have been made to support Mac OS X 10.2. Shared
libraries are not supported on this architecture, but everything else
should work.
(WB 2002/12/18)
</p>
- <li> <p>
+ <li> <p>
New: deal.II can be compiled with version 7.0 of
Intel's icc compiler, which was recently released. Since this compiler
finally supports the very restrictive flags <code>-Xc -ansi</code> that
(WB 2002/12/05)
</p>
- <li> <p>
+ <li> <p>
Changed: Previously, we just set the preprocessor variable
<code>DEAL_II_USE_MT</code>, when <code>--with-multithreading</code> was
given as argument to <code>./configure</code>. Tests in the code
therefore looked like <code>#ifdef DEAL_II_USE_MT</code>. This has been
changed so that the variable is always defined, but its value is now
- equal to <code>1</code>
+ equal to <code>1</code>
when multithreading was requested, and zero otherwise. The reason for
this is that you can now write <code>if (DEAL_II_USE_MT && ...)</code>
conditions, and need not interleave if-else clauses from regular code
(WB 2002/11/14)
</p>
- <li> <p>
+ <li> <p>
New: Object files are now named according to the local defaults
on the system we are running on. On Unix systems, this is
usually a <code>.o</code> suffix, while on Windows it is
(WB 2002/11/11)
</p>
- <li> <p>
+ <li> <p>
New: deal.II can now also be compiled with Red Hat's
version of the gcc compiler, gcc 2.96. However, some problems remain;
for more information.
(WB 2002/10/14)
</p>
- <li> <p>
+ <li> <p>
Fixed: On CygWin, one header files has a <code>#define quad
quad_t</code>. This is annoying, since we have local variables and
member functions with the name <code>quad</code>, and in fact it breaks
set.
<br>
(Stephen Kolaroff, WB
- 2002/09/28)
+ 2002/09/28)
</p>
- <li> <p>
+ <li> <p>
New: Since <code>math.h</code> only defines the values of PI or E (as
<code>M_PI</code> and <code>M_E</code>) when certain defines are set (on
Linux, these are <code>__USE_BSD</code> or <code>__USE_XOPEN</code>),
(WB 2002/09/12)
</p>
- <li> <p>
+ <li> <p>
New: <code>base/config.h</code> now exports the
deal.II base directory through the
<code>DEAL_II_PATH</code> preprocessor variable.
(WB 2002/09/09)
</p>
- <li> <p>
+ <li> <p>
Removed: The forward declarations files have gone. We have never
propagated their use in the example programs, but these files have been
there in the base, lac, and grid include directories, and forward
(WB 2002/09/05)
</p>
- <li> <p>
+ <li> <p>
New: There is now
a new report
- on assembling matrices available from the
+ on assembling matrices available from the
<a href="../documentation.html" target="body">documentation
page</a>. The main focus is assembling of matrices for
vector-valued problems, where shape functions are
vector-valued, and may have only one or more non-zero vector
- components.
+ components.
<br>
(WB 2002/06/17)
</p>
(WB 2002/06/10)
</p>
- <li> <p>
+ <li> <p>
New: deal.II now also supports vector-valued
finite elements with shape functions for which more than just
one vector component is non-zero. Such elements are, for
(WB 2002/06/10)
</p>
- <li> <p>
+ <li> <p>
New: The top-level makefile now how a target <tt>distclean</tt>.
<tt>clean</tt> leaves the libraries now, removing
everything that is not needed to use
(WB 2003/02/13)
</p>
- <li> <p>
+ <li> <p>
New: The <code>Function</code> class now has an assignment
operator. This way, you can put function objects into
containers. However, the assignment operator makes sure that only
(WB 2003/02/06)
</p>
- <li> <p>
+ <li> <p>
New: The <code>ThreadMutex</code> classes now have a
member class <code>ScopedLock</code> that implements the
scoped thread-safe locking pattern of Doug Schmidt. It is also used in
(WB 2003/01/28)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code>PosixThreadManager</code> called
its <code>wait</code> function in the
destructor. If this had been called before already, then the
(Michael Anderson, WB 2003/01/16)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code>Subscriptor</code> uses a counter to
count how many <code>SmartPointer</code> objects subscribe
to the pointed-to object. This counter needs to be a volatile variable
(WB 2003/01/11)
</p>
- <li> <p>
+ <li> <p>
Fixed: In multithreaded mode, when a new thread is started, the
arguments to the function being called need to be copied from the stack
of the starting thread to that of the new thread. In order to
(Michael Anderson, WB 2003/01/11)
</p>
- <li> <p>
+ <li> <p>
New: There are now classes <code>ThreadCondition</code>
that implement thread condition variable operations through POSIX
functions (in multithreaded mode) or doing nothing (in singlethreaded
(WB 2003/01/11)
</p>
- <li> <p>
+ <li> <p>
New: Newer versions of gcc have a very nice feature: you can set
a verbose terminate handler, that not only aborts a program
when an exception is thrown and not caught somewhere, but
(WB 2002/12/19)
</p>
- <li> <p>
+ <li> <p>
New: There is now a <code>Patterns::List</code> pattern
that matches a list of elements each of which has to satisfy a pattern
that was given to the constructor.
(WB 2002/11/29)
</p>
- <li> <p>
+ <li> <p>
Changed: In POSIX mode, when the <code
class="member">ThreadManager</code> class created a new thread through
<code>pthread_create</code>, it only checked for the
(WB 2002/11/13)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code>write_text</code> and <code
class="member">write_tex</code> functions of the <code
class="class">TableHandler</code> class now check whether their
(RH 2002/11/11)
</p>
- <li> <p>
+ <li> <p>
New: Added Hierarchical Polynomial (similar to Legendre class). Will
eventually be used in a hierarchical FiniteElement class similar to
FE_Q class. Included in Polynomials namespace.
(RH 2002/09/24)
</p>
- <li> <p>
+ <li> <p>
New: The <code>Threads::n_existing_threads</code>
function returns the present number of existing threads, allowing an
assessment whether it is useful to spawn new threads or rather perform
(GK 2002/09/08)
</p>
- <li> <p>
+ <li> <p>
Extended: Previously, the <code
class="class">Threads::PosixThreadBarrier</code>
class could not be used at all (they threw exceptions), if your system
(WB 2002/09/05)
</p>
- <li> <p>
+ <li> <p>
New: The old class <code>vector2d</code>, implementing a
two-dimensional array of objects is now gone. Instead, there is the new
<code>Table</code> class that implements tables of
(WB 2002/09/03)
</p>
- <li> <p>
+ <li> <p>
New: There are now functions returning the transpose of <code
class="class">Tensor</code> objects of rank 2.
<br>
(WB 2002/08/31)
</p>
- <li> <p>
+ <li> <p>
New: Row accessors for the <code>vector2d</code>
class now have a member function <code
class="member">size</code> that returns the size of the row,
(WB 2002/07/24)
</p>
- <li> <p>
+ <li> <p>
Fixed: In EPS output, colors were set to invalid values if the
values of the field that is used for coloring are all
equal. This happens, for example, in the very first time step
(WB 2002/07/24)
</p>
- <li> <p>
+ <li> <p>
Changed: To save disk space, color values in EPS output are
written as grey scale with only one value instead of three RGB
values if the color so represented is actually a grey scale.
(WB 2002/07/24)
</p>
- <li> <p>
+ <li> <p>
New: There are now operators to allow multiplication and
division of <code>Tensor</code> objects by scalar
factors.
<h3>lac</h3>
<ol>
- <li> <p> New: Function
+ <li> <p> New: Function
<code>BlockSparseMatrix::print_formatted</code>
- added, similar to the
+ added, similar to the
<code>BlockVector::print_formatted</code> function.
<br>
(Brian Carnes 2003/06/02)
- </p>
+ </p>
<li> <p> New: Functions <code>SparseMatrix::operator *=</code>
- and <code>SparseMatrix::operator /=</code> as well as
+ and <code>SparseMatrix::operator /=</code> as well as
<code>BlockSparseMatrix::operator *=</code>
- and <code>BlockSparseMatrix::operator /=</code>
+ and <code>BlockSparseMatrix::operator /=</code>
are added.
<br>
(Brian Carnes 2003/06/02)
- </p>
+ </p>
<li> <p> Deprecated: The functions <code>Vector::scale</code>
and <code>BlockVector::scale</code> are now deprecated
<code>operator/=</code> instead.
<br>
(WB 2003/05/31)
- </p>
+ </p>
<li> <p> New: <code>Vector</code>, <code
class="class">BlockVector</code> and
for scaling by a scalar.
<br>
(WB 2003/05/31)
- </p>
+ </p>
<li> <p> New: <code>PointerMatrix</code> now
has <code>empty()</code> function, which returns true if the pointer is
(WB 2003/02/13)
</p>
- <li> <p>
+ <li> <p>
New: The <code>SolverGMRES</code> now accepts an
<code>AdditionalData</code> parameter
<code>use_default_residual</code> whose default is
(RH 2003/01/31)
</p>
- <li> <p>
+ <li> <p>
New: <code>FullMatrix</code> has a function <code
class="member">copy_from</code>, copying from sparse matrices.
It uses iterators of the sparse matrix classes.
<br>
(Stephen Kolaroff 2002/11/27)
</p>
-
+
<li> <p>
Changed: In multithread mode, the <code
class="class">SparseMatrix</code> would spawn
<br>
(WB 2002/11/13)
</p>
-
+
<li> <p>
- Fixed: In the <code>SparseDirectMA27</code> class, wrapping
+ Fixed: In the <code>SparseDirectMA27</code> class, wrapping
the MA27 solver written in Fortran77 into some structure amenable to C++,
we wrote beyond the end of an array on rare inputs. This is now fixed. The
same probably holds for the respective code for MA47.
computations using this solver. The actual operations are therefore
distributed to distint programs that have separate address spaces. This
allows to have as many concurrent instances of this solver in parallel
- as you want. For more information, read the documentation of the
+ as you want. For more information, read the documentation of the
<code>SparseDirectMA27</code> class.
<br>
(WB 2002/09/25)
<li> <p>
New: Classes <code>SparsityPattern</code> and <code
- class="class">SparseMatrix</code> now have functions <code
+ class="class">SparseMatrix</code> now have functions <code
class="member">block_write/block_read</code>, allowing to dump the data
of these objects into a file in binary format, and later to re-read it
without much need for parsing.
(GK 2002/08/13)
</p>
- <li> <p>
+ <li> <p>
New: a way of using abstract base classes for matrices has
been implemented with <code>PointerMatrixBase</code>
and <code>PointerMatrix</code>. Storing a matrix in
to the child cell; <code
class="member">child_to_cell_coordinates</code> that does
exactly the opposite; and <code
- class="member">is_inside_unit_cell</code> that tells whether a
+ class="member">is_inside_unit_cell</code> that tells whether a
given point is inside the unit cell.
<br>
(WB 2003/05/01)
(WB 2003/02/21)
</p>
- <li> <p>
+ <li> <p>
Changed: Some implementation details of the <code
class="class">GridReordering</code> class have been moved to a
namespace <code>internal</code>.
</p>
<li> <p>
- Fixed: Objects of type <code>FESystem</code>
+ Fixed: Objects of type <code>FESystem</code>
could not be constructed if one of the elements it is to be composed of
did not implement interface constraints for hanging nodes. This is now
fixed: you can construct such a composed element, but it does not
(RH 2003/02/03)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code>GridReordering</code> tried to be
thread-safe in the initialization of some data, but was not due to a
typo. This is now fixed.
(WB 2003/01/28)
</p>
- <li> <p>
+ <li> <p>
Changed: The <code>FEValues::get_cell</code> and
<code>FEValues::get_face</code> functions have
been removed, since they limited our ability to use this
(WB 2003/01/17)
</p>
- <li> <p>
+ <li> <p>
New: The DoF accessor classes now have a function <code
class="member">get_fe()</code> that returns a reference to the finite
element being used on this cell. The result is of course identical to
(WB 2003/01/10)
</p>
- <li> <p>
- New: Checked in new <code>GridGenerator</code>
+ <li> <p>
+ New: Checked in new <code>GridGenerator</code>
member function <code>half_hyper_ball</code>,
- derived from member <code>hyper_ball</code>.
+ derived from member <code>hyper_ball</code>.
The initial mesh contains four elements. This mesh will work with
the boundary class <code>HalfHyperBallBoundary</code>.
<br>
(Brian Carnes 2002/12/16)
</p>
- <li> <p>
+ <li> <p>
New: Checked in new class <code>FE_Q_Hierarchical</code>
- derived from class <code>FiniteElement</code>.
- This element is analogous to <code>FE_Q</code>, but
- makes use of hierarchical shape functions, based on the
- <code>Polynomials::Hierarchical</code> class.
+ derived from class <code>FiniteElement</code>.
+ This element is analogous to <code>FE_Q</code>, but
+ makes use of hierarchical shape functions, based on the
+ <code>Polynomials::Hierarchical</code> class.
For <code>degree>1</code>, the non-nodal basis functions are "bubble"
functions, which are not Lagrange polynomials. Therefore, the usual
interpolation based on using unit support points will not work for
interpolation-projection operator, based on an hp-type interpolant.
<br>
The files for this element are
- <code>deal.II/include/fe/fe_q_hierarchical.h</code> and
+ <code>deal.II/include/fe/fe_q_hierarchical.h</code> and
<code>deal.II/source/fe/fe_q_hierarchical.cc</code>.
<br>
(Brian Carnes 2002/12/13)
</p>
- <li> <p>
+ <li> <p>
New: For finite element classes, the functions
<code>unit_support_point</code> and
<code>unit_face_support_point</code> return the position
of an individual support point. This is necessary when you want to get
information about the support points of certain components in a composed
finite element, where not all components provide support points, and the
- composed element thus does not fill the array the
+ composed element thus does not fill the array the
<code>get_unit_support_points</code> function returns.
<br>
(WB 2002/11/23)
</p>
- <li> <p>
+ <li> <p>
Fixed: Vectors could not be given as input and output vectors to the
<code>SolutionTransfer</code> class at the same time, but
this was not checked. An assertion has now been added to ensure this
- requirement.
+ requirement.
<br>
(WB 2002/10/28)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code
class="member">DoFRenumbering::component_wise</code> function accepts a
parameter indicating the order in which the degrees of freedom
(WB 2002/10/17)
</p>
- <li> <p>
+ <li> <p>
New: The <code
class="member">FiniteElement::system_to_base_index</code> function now
exports the values of the <code
(WB 2002/10/17)
</p>
- <li> <p>
+ <li> <p>
New: The <code
class="member">FiniteElement::element_multiplicity</code> function was
previously only available in the <code>FESystem</code>
(WB 2002/10/17)
</p>
- <li> <p>
+ <li> <p>
Fixed: The algorithm to generate neighbor information for patches from
cells in <code>DataOut::build_patches</code> was rather
inefficient, tripling the time for patch generation when support for
(WB 2002/10/10)
</p>
- <li> <p>
+ <li> <p>
Changed (internals): Previously, the finite element base class
initialized the restriction, prolongation, and face constraints matrices
to their correct size. Derived classes had to fill these classes, and
(GK 2002/09/19)
</p>
- <li> <p>
+ <li> <p>
Fixed: Some of the various instances of the <code
class="member">VectorTools::interpolate_boundary_values</code> functions
were not explicitly instantiated, leading to linker errors. This is now
- fixed.
+ fixed.
<br>
(WB 2002/09/19)
</p>
- <li> <p>
+ <li> <p>
Removed: The <code
class="member">FiniteElement::component_to_system_index</code> function
and its counterpart for faces is gone. This data did not make much sense
(WB 2002/09/16)
</p>
- <li> <p>
+ <li> <p>
Changed: The <code
class="member">FiniteElement::restriction_is_additive</code> function
used to take an argument that denoted the vector component of a finite
(WB 2002/09/16)
</p>
- <li> <p>
+ <li> <p>
New: The mapping classes now also know how to transform tensors of rank
2, not only vectors (tensors of rank 1) in a co- and contravariant way.
<br>
(WB 2002/09/03)
</p>
- <li> <p>
+ <li> <p>
Fixed: the <code>GridIn</code> class had problems
when reading in UCD grids with comment lines that contained
only the comment sign, but nothing else. This is now fixed.
(WB 2002/08/30)
</p>
- <li> <p>
+ <li> <p>
Improved: <code>VectorTools</code>::<code
class="member">integrate_difference</code> can compute <i>L<sup>p</sup></i>
and <i>W<sup>1,p</sup></i> norms for arbitrary <i>p</i>. The function
the library.
<br>
(RH 2002/06/14)
- </p>
+ </p>
<li> <p> Fixed: <code
class="member">FE_Q::has_support_on_face</code> always returned
(RH 2002/06/10)
</p>
- <li> <p>
+ <li> <p>
Fixed: The <code>get_dof_indices</code>
functions of DoF accessor classes used to work also for
non-active cells. However, the results were bogus except for
(WB 2002/06/08)
</p>
- <li> <p>
+ <li> <p>
New: For encapsulated postscript output of 2d grids, it is now
possible to tell the <code>GridOut</code> class to
write the cell numbers into each cell, as well as the numbers
/**
* @page changes_between_4_0_and_5_0 Changes between Version 4.0 and 5.0
-
+
<p>
This is the list of changes made between the deal.II releases listed above.
All entries are signed with the names of the author. Regular
that take arguments. Instead, you should now use <code>matrix=0</code> to simply set all elements of the
matrix to zero. If you want to notify a sparse matrix that its sparsity
pattern has changed, use the <code>reinit(SparsityPattern)</code> function.
- <br>
+ <br>
(WB 2004/05/10)
</p>
recommended way is to use the more obvious notation <code>vector=0</code> and <code>matrix=0</code>. To
reset the elements of a table over arbitrary objects, use
<code>Table<T>::reset_values()</code>.
- <br>
+ <br>
(WB 2004/05/10)
</p>
<li> <p> Removed: The <code>SparseLUDecomposition::reinit</code> and <code>SparseMIC::reinit</code> functions without
argument had been deprecated before, and have now been removed.
- <br>
+ <br>
(WB 2004/05/10)
</p>
therefore been renamed <code
class="member">get_restriction_matrix</code>, which also better
reflects what it is actually doing. To keep consistent, we have also
- rename the corresponding function <code>prolongate</code> to <code>get_prolongation_matrix</code>.
+ rename the corresponding function <code>prolongate</code> to <code>get_prolongation_matrix</code>.
<br>
(WB 2004/02/29)
</p>
these quadrature formula onto the faces from the other side. We need
this due to the fact that we now support faces in 3d that have a normal
vector opposite to the standard direction.
- <br>
+ <br>
(WB 2003/10/17)
</p>
<li> <p> Moved and changed: The header file
- <tt>include/numerics/dof_renumbering.h</tt> has been moved to the
+ <tt>include/numerics/dof_renumbering.h</tt> has been moved to the
directory <tt>include/dofs</tt>, where it logically
belongs. Furthermore, the sorting parameter of the function <code>DoFRenumbering</code><code>::component_wise</code> has changed its meaning. See
the reference documentation for details.
(WB 2004/04/12)
</p>
- <li> <p> New: deal.II is now able to interface to the
+ <li> <p> New: deal.II is now able to interface to the
<a href="http://www-users.cs.umn.edu/~karypis/metis/index.html"
target="_top">METIS</a> library to generate domain partitionings. This
is enabled if a METIS installation is detected, which either happens
automatically in <code>./configure</code>, or
by passing a value to the switch <code>--with-metis</code> to configure
the path of a METIS installation. For more information see the README
- file.
+ file.
<br>
(WB 2004/03/08)
</p>
<li> <p> New: deal.II now comes with a complete set of
wrappers classes for <a href="http://www.mcs.anl.gov/petsc/"
- target="_top">PETSc</a> vectors, matrices, linear solvers and
+ target="_top">PETSc</a> vectors, matrices, linear solvers and
preconditioners. Many of the algorithms in deal.II have also been
updated to make use of these wrappers. All of this is only enabled if a
PETSc installation is detected. This either happens automatically in
some compiler/platforms requires boost files in addition to
those included in the <tt>contrib</tt> directory please report
the missing boost files to the deal.II mailing list.
- <br>
+ <br>
(RH 2004/02/16)
</p>
Changed: We don't support compilation by Intel'c icc compiler version 5
anymore. Since the present version of this compiler is 8, this is
probably not a real problem.
- <br>
+ <br>
(WB 2003/12/20)
</p>
Fixed: <code>step-9</code> had the computation of the value of the
solution in the mid-point of the neighbor cell wrong. This is now
fixed. Consequently, the resulting mesh also looks much nicer now (much
- smoother).
- <br>
+ smoother).
+ <br>
(<a href="mailto:werner.scholz@tuwien.ac.at">Werner Scholz</a>
2003/12/11)
</p>
unsigned integer. This value is widely used throughout the library as a
marker for an invalid unsigned integer value, such as an invalid array
index, an invalid array size, and the like.
- <br>
+ <br>
(WB 2003/10/24)
</p>
<li> <p>
Augmented: The <code
- class="member">GeometryInfo::child_cell_on_face</code>
+ class="member">GeometryInfo::child_cell_on_face</code>
results in a result that might not be what you expect in 3d in some
cases. A warning has been added to the documentation, as well as a
reference to the function to call instead.
- <br>
+ <br>
(WB 2003/10/17)
</p>
(WB 2003/09/23)
</p>
- <li> <p> Fixed: On AIX systems, the xlf Fortran77 compiler wasn't recognized
+ <li> <p> Fixed: On AIX systems, the xlf Fortran77 compiler wasn't recognized
due to the fact that when called with -v it generates many pages
of output, later leading to a shell error. This is now fixed.
<br>
<h3>base</h3>
<ol>
- <li> <p>
+ <li> <p>
New: There is now a new <code>PolynomialsP</code>
class which is based on <code>Polynomials::Monomial</code> and <code>PolynomialSpace</code>. In contrast to the
default ordering of the polynomials in <code>PolynomialSpace</code>, (e.g. for degree=2) <i>1,
(RH 2004/03/11)
</p>
- <li> <p>
+ <li> <p>
New: The classes <code>PolynomialSpace</code> and
<code>TensorProductPolynomials</code> now have
new <code>set_numbering</code> functions which
<li> <p> Fixed: The constructor of the <code>QGauss</code> class
computed positions and weights of quadrature points in long double accuracy.
- However, on machines where long double is the same as double, it
+ However, on machines where long double is the same as double, it
never reached the requested accuracy, in effect leading to an infinite loop.
This is now fixed.
<br>
<br>
(WB 2004/01/26)
</p>
-
+
<li> <p> Improved: Initialization routines of class <code>SparseMatrix</code> have an additional parameter
controlling the storage of diagonal entries.
<br>
(GK 2003/11/18)
</p>
- <li> <p> New:
+ <li> <p> New:
<code>SolverFGMRES</code> implements the flexible
GMRES method with varying preconditioner from the right. It is
also accessible in <code>SolverSelector</code> by choosing <tt>fgmres</tt>.
<br>
(Luca Heltai 2004/04/21)
</p>
-
+
<li> <p>
New: The function <code>GridGenerator::cylinder_shell</code> generates a domain
of the type which the name suggests.
<br>
(WB 2004/04/19)
</p>
-
+
<li> <p>
Changed: The <code>KellyErrorEstimator::estimate</code> function takes an
additional parameter that lets it only compute error indicators for a
</p>
<li> <p> New: The new functions <code>SparsityPattern::partition</code>, <code>GridTools::partition_triangulation</code>, <code>DoFTools::get_subdomain_association</code>, <code>DoFTools::count_dofs_with_subdomain_association</code>,
- <code>GridTools::get_subdomain_association</code>, <code>GridTools::count_cells_with_subdomain_association</code>,
+ <code>GridTools::get_subdomain_association</code>, <code>GridTools::count_cells_with_subdomain_association</code>,
and <code>DoFRenumbering::subdomain_wise</code> can now
be used to generate partitions of a triangulation and its associated
degrees of freedom suitable for parallel computations with PETSc.
Improved: <code>FiniteElementData</code> has a function
<code>tensor_degree()</code>, returning the degree of the
polynomial space suitable for choosing a tensor product quadrature
- formula.
+ formula.
<br>
(GK 2003/11/28)
</p>
<li> <p>
New: There are now functions <code>CellAccessor::neighbor_child_on_subface</code> and <code>DoFCellAccessor::neighbor_child_on_subface</code>
that should be called instead of using <code>GeometryInfo::child_cell_on_face</code> in most cases.
- <br>
+ <br>
(WB 2003/10/17)
</p>
<li> <p>
New: <code>GridGenerator</code> has a new
- function <code>subdivided_hyper_rectangle</code>
- which generates a rectangle with given corner points and possibly
+ function <code>subdivided_hyper_rectangle</code>
+ which generates a rectangle with given corner points and possibly
different numbers of subdivisions in different directions.
Use it, e.g., to generate a domain of 1*4 length units
with square cells.
- <br>
+ <br>
(Joerg Weimar 2003/09/09)
</p>
Improved: The 3d grid reordering code in the <code>GridReordering</code> class now uses an algorithm
that is linear in the number of elements. The old code was
exponential, so this is a vast improvement.
- <br>
+ <br>
(Michael Anderson 2003/09/23)
</p>
<li> <p>
Improved: <code>GridOut</code> has a an improved
- functionality for <code>write_eps</code>
+ functionality for <code>write_eps</code>
to color the grid according to the refinement level.
- A corresponding option is included in
+ A corresponding option is included in
<code>GridOutFlags::Eps<2></code>.
- <br>
+ <br>
(Joerg Weimar 2003/09/09)
</p>
Improved: <code>FEFaceValuesBase</code> has a new
function <code>orientation</code> accessing a unique
and consistent orientation for each face.
- <br>
+ <br>
(GK 2003/06/23)
</p>
- <li> <p>
+ <li> <p>
Changed: Embedding and restriction matrices for intergrid transfer are
now computed in the constructor of most elements, rather than taken from
precomputed and tabulated values. This removes restrictions on which
(WB 2003/06/09)
</p>
- <li> <p>
+ <li> <p>
New: Finite elements got a function <code>get_interpolation_matrix</code>, with which they can
compute interpolations between different finite elements. Most will use
this to compute interpolations between finite elements of the same kind
(WB 2003/06/09)
</p>
- <li> <p>
+ <li> <p>
New: Finite elements got a function <code>get_name</code>, which can be used to identify a finite
element by its name.
<br>
(WB 2003/06/09)
</p>
- <li> <p>
+ <li> <p>
New: Raviart-Thomas elements are now implemented in the <code>FE_RaviartThomas</code> class.
<br>
(WB 2003/06/09)
and leads to very hard to find errors. Therefore, all constructors and
reinit functions of the parallel sparse matrix classes now take an
additional argument indicating the partitioning of columns of the matrix.
- <br>
+ <br>
(WB 2004/06/02)
</p>
</ol>
<ol>
<li> <p>
- New: Support of shared libraries under Cygwin / Windows systems. Shared
+ New: Support of shared libraries under Cygwin / Windows systems. Shared
libraries are used as default. To use, make sure the library path is included
in your <code>PATH</code> environment variable.
<br>
</p>
<li> <p>
- New: Configuration now detects the Intel Fortran compiler and can set
+ New: Configuration now detects the Intel Fortran compiler and can set
compilation flags accordingly.
<br>
(WB 2004/11/04)
</p>
<li> <p>
- New: configuration option <tt>--with-umfpack</tt> for using the
+ New: configuration option <tt>--with-umfpack</tt> for using the
UMFPack version and enabling the class <code
class="class">SparseDirectUMFPACK</code>, both included by
Wolfgang.
<li> <p>
Fixed: None of the formulas in the step-8 tutorial program web
page were visible. This is now fixed.
- <br>
+ <br>
(WB 2004/06/29)
</p>
<h3>base</h3>
<ol>
- <li> <p>
+ <li> <p>
New: The <code>MultithreadInfo</code> class now also detects multiple
processors on Mac OS X.
- <br>
+ <br>
(Helmut Müller 2004/11/29)
</p>
- <li> <p>
+ <li> <p>
New: The <code>
TableHandler::write_tex</code> now accepts the additional boolean
argument <code>with_header</code> which is set to
- true by default and tells the function whether to add the latex
+ true by default and tells the function whether to add the latex
header and footer (i.e. the \\documentclass{...},
- \\begin{document} and \\end{document} stuff) to the table.<br>
- In addition to this, there are two new members in the above class:
+ \\begin{document} and \\end{document} stuff) to the table.<br>
+ In addition to this, there are two new members in the above class:
<code> TableHandler::tex_set_table_caption</code> and
<code> TableHandler::tex_set_table_label</code> to
- add a caption and a label to the tex generated table.
- <br>
+ add a caption and a label to the tex generated table.
+ <br>
(Luca Heltai 2004/10/29)
</p>
class="class">Table<N,T></code> class represents an
<code>N</code>-dimensional array and might replace the
N-times-nested-use of the <code
- class="class">std::vector</code> class.
+ class="class">std::vector</code> class.
<br>
(RH 2004/08/13)
</p>
<li> <p>
New: Class <code>TableIndices</code> now has
operators that check for equality and inequality of objects.
- <br>
+ <br>
(WB 2004/07/28)
</p>
<li> <p>
New: A class <code>PointerComparison</code> for comparing
pointers that may or may not be of the same type.
- <br>
+ <br>
(WB 2004/06/22)
</p>
class="member">ok</code> function. It was never really used for
anything, and has thus finally been removed. The <code
class="member">ok</code> is consequently gone as well.
- <br>
+ <br>
(WB 2004/06/06)
</p>
or not. This is particular useful for parallel computations
when only one process should print to standard output, while in all
other processes we simply want output suppressed.
- <br>
+ <br>
(RH 2004/05/26)
</p>
</ol>
prevented deal.II from working with it (it couldn't be compiled at
all). We have added checks that make sure that it can now again be
compiled with PETSc versions 2.1.6, 2.2.0, and 2.2.1.
- <br>
+ <br>
(WB 2004/10/07)
</p>
<li> <p>
New: A class <code>PreconditionLU</code> which
provides a wrapper to the complete LU decomposition
- preconditioner of PETSc. Furthermore a class
+ preconditioner of PETSc. Furthermore a class
<code>SolverPreOnly</code> was implemented. It
is a wrapper for the PETSc solver type KSPPREONLY, which
only applies the preconditioner. In conjunction with
<code>PreconditionLU</code> this provides a
simple direct solver, which could be used for small to medium
sized problems on a single processor machine.
- <br>
+ <br>
(Oliver Kayser-Herold 2004/07/27)
</p>
-
+
<li> <p>
Improved: <code>VectorTools::point_difference</code>
used to use an algorithm to find the cell the given point is in that
to use
<code>GridTools::find_active_cell_around_point</code>
which is only logarithmic in its complexity.
- <br>
+ <br>
(WB 2004/07/07)
</p>
(and invalid) states when some of the blocks had empty rows. In this
common case, we would frequently skip elements when looping over the
elements of a block matrix. These cases should now be fixed.
- <br>
+ <br>
(WB 2004/07/07)
</p>
the position within a row. However, this was fragile, and has been
removed. If you want an ordering of elements within a row, use the
<code>operator <</code> to compare iterators.
- <br>
+ <br>
(WB 2004/07/07)
</p>
class="class">PETScWrappers::SolverGMRES::AdditionalData</code>
class now takes an additional flag indicating the use of left
or right preconditioning.
- <br>
+ <br>
(RH 2004/06/24)
</p>
<li> <p>
Improved: The <code>SparseDirectMA27</code> class can now
handle float as well as double input matrices and vectors.
- <br>
+ <br>
(WB 2004/06/23)
</p>
and output vectors, but were only instantiated in case the arguments
were the same, and could also not be compiled if they weren't. Both
problems are now fixed.
- <br>
+ <br>
(WB 2004/06/22)
</p>
average. This should significantly reduce the total amount of memory
required as well as memory fragmentation. It also cuts run-time for
element addition by more than half.
- <br>
+ <br>
(WB 2004/06/21)
</p>
class="member">CompressedSparsityPattern::max_entries_per_row()</code>
ignored the first row and thus sometimes returned a value that was too
low. This is now fixed.
- <br>
+ <br>
(WB 2004/06/21)
</p>
removed them, but it still kept the number of blocks
unchanged. Subsequent accesses to these blocks, or to their information
yielded segmentation faults.
- <br>
+ <br>
(WB 2004/06/02)
</p>
still exist in all their functionality, this allowed us to now also
have block matrices and vectors for PETSc sequential and parallel
objects.
- <br>
+ <br>
(WB 2004/06/02)
</p>
temporaries without our ado. This bug is now detected during
configuration time of the library, and these constructors are available
whenever the compiler does not contain this bug.
- <br>
+ <br>
(WB 2004/06/02)
</p>
now also have member functions <code>local_range</code>,
<code>in_local_range</code>, and <code
class="member">local_size</code>.
- <br>
+ <br>
(WB 2004/06/02)
</p>
New: Parallel PETSc matrix and vector classes now have member functions
<code>get_mpi_communicator</code> that returns the MPI
communicator object these objects operate on.
- <br>
+ <br>
(WB 2004/06/02)
</p>
than a nonconstant reference to the MPI communicator to be used. This
prevents some unnecessary compiler problems in conjunction with the new
<code>get_mpi_communicator</code> function.
- <br>
+ <br>
(WB 2004/06/02)
</p>
initialized (via either constructor or reinit functions) with a
compressed sparsity pattern object, allowing for more efficient
preallocation of nonzero entries.
- <br>
+ <br>
(WB 2004/06/02)
</p>
</ol>
places in the library use or touch the user pointers. That is
wrong, the <code>SolutionTransfer</code> class
actually does. This is now properly documented.
- <br>
+ <br>
(WB 2004/09/15)
</p>
class="member">DerivativeApproximation::approximate_second_derivative</code>
function produced wrong results ("not a number", or an exception instead of
"zero") if the field it was given was constant. This is now fixed.
- <br>
+ <br>
(WB 2004/08/05)
</p>
New: The <code>MatrixTools::apply_boundary_values</code>
that works on block matrices and vectors is now templatized over the
number type, i.e. it also works for float matrices and vectors.
- <br>
+ <br>
(WB 2004/06/22)
</p>
New: The new <code>FEValuesBase::get_cell</code>
function returns the current cell, i.e. the latest cell the
<code>FEValues</code> object was reinited with.
- <br>
+ <br>
(RH 2004/06/22)
</p>
of information, but the cases where this applies and the strategies to
work with this problem anyway are now well documented in the
documentation of the <code>MatrixTools</code> class.
- <br>
+ <br>
(WB 2004/06/07)
</p>
Jacobians. In usual finite element computations this simply leads to
all components of the linear system being negated, so it is not
harmful. It is fixed now anyway.
- <br>
+ <br>
(WB 2004/06/07)
</p>
</ol>
<li> <p> New:
The step-19 tutorial demonstrates handling of parameters from a
parameter file, as well as some simple techniques for merging output
- from parallel computations.
+ from parallel computations.
<br>
(WB 2005/09/09)
</p>
might issue error messages about multiple defined symbols. This is now
detected, and the linker is forced to ignore these errors. However, if
you accidentally defined symbols twice you might have a hard time
- debugging now... In this case remove the <code>-Xlinker
+ debugging now... In this case remove the <code>-Xlinker
--allow-multiple-definition</code> flag from Make.global_options.
<br>
(Ralf B. Schulz, 2005/07/13)
<li> <p>
Fixed: When compiling shared libraries on CygWin systems, warnings
concerning the <code>-fPIC</code> option were issued by the compiler.
- This is now fixed.
+ This is now fixed.
Also, configure now issues a message at the end that you should include
the DLL-path in your <code>.bash_profile</code> file on these systems.
<br>
systems. A new preprocessor variable, <code>DEAL_II_BROKEN_SOCKETS</code>
has been added to <code>base/config.h</code> which is defined on
affected systems.
- <br>
+ <br>
(Ralf B. Schulz, WB, 2005/03/02)
</p>
Fixed: The step-16 example program wasn't listed in the
navigation bar of the tutorial section, although it was in the
table of contents. This is fixed now.
- <br>
+ <br>
(WB, 2005/02/09)
</p>
</ol>
a file in a sequence of directories and by appending suffixes. The class
generates the complete file name including directory prefix and suffix
such that it can be used for subsequently opening the file.
- <br>
+ <br>
(GK, 2005/09/14)
</p>
<br>
(Ralf B. Schulz, 2005/07/19)
</p>
-
+
<li> <p>
Removed: The <code>write_multigrid</code> flag in <code
- class="member">DataOutBase::DXFlags</code>
+ class="member">DataOutBase::DXFlags</code>
has been removed, since it wasn't used anywhere.
<br>
(WB, 2005/07/14)
<br>
(Ralf B. Schulz, 2005/07/13)
</p>
-
+
<li> <p>
Improved: The <code>QProjector</code> now has
functions <code>project_to_face</code> and <code
class="member">double_contract</code> contracts two tensors of
rank 4 and 2 into a tensor of rank 2, by contracting over two
indices at the same time.
- <br>
+ <br>
(WB, 2005/03/29)
</p>
class="member">TableIndicesBase::sort</code> allows to sort the indices
in ascending order. This is useful for cases where the order of indices
is irrelevant (for example in symmetric tables).
- <br>
+ <br>
(WB, 2005/03/29)
</p>
<li> <p>
New: There is a new class <code>SymmetricTensor</code>
that provides storage and operations for symmetric tensors.
- <br>
+ <br>
(WB, 2005/03/28)
</p>
easily. The additional feature is being incorporated into <code
class="class">SmartPointer</code> constructors throughout the
library.
- <br>
+ <br>
(GK, 2005/03/16)
</p>
<li> <p>
New: Class <code>FunctionParser</code>. Wrapper
- class for the fparser library (see
+ class for the fparser library (see
<a href="http://warp.povusers.org/FunctionParser/">
http://warp.povusers.org/FunctionParser/</a>).
- <br>
+ <br>
(Luca Heltai, 2005/03/07).
</p>
<li> <p>
Fixed: The class <code
- class="class">MultipleParameterLoop::UserClass</code> had only
+ class="class">MultipleParameterLoop::UserClass</code> had only
virtual abstract functions but no virtual destructor. This caused
warnings with some compilers, and is generally bad practice
anyway. This is now fixed. The same holds with respect to the class
<code>DataOutInterface</code>.
- <br>
+ <br>
(WB, 2005/02/20)
</p>
</ol>
New: The new function <code
class="member">PetscWrappers::VectorBase::set</code> allows to set
several vector elements at once.
- <br>
+ <br>
(WB, 2005/08/10)
</p>
class="member">clear_rows</code> that allow to set the elements of a row
of a matrix to zero, without having to traverse the individual
elements.
- <br>
+ <br>
(WB, 2005/08/10)
</p>
has the complete interface of a preconditioner
class. Therefore, it can be used in iterative solvers and in
multigrid smoothers.
- <br>
+ <br>
(GK, 2005/05/10)
</p>
<li> <p>
Fixed: The PETSc matrix iterators had trouble when some rows of a
matrix were empty. This has now been mostly fixed.
- <br>
+ <br>
(WB, 2005/05/02)
</p>
about its effectiveness. This has now been fixed: The code now properly
initializes these elements, and makes the resulting matrix much faster
to use.
- <br>
+ <br>
(WB, 2005/04/04)
</p>
New: The <code>ProductSparseMatrix</code>
implements the product of two rectangular sparse matrices with
the same <code>value_type</code>
- <br>
+ <br>
(GK, 2005/03/11)
</p>
<li> <p>
New: The <code>PreconditionRichardson</code>
implements a Richardson preconditioner.
- <br>
+ <br>
(GK, 2005/03/10)
</p>
made it a little awkward to use in some places. This has now
been fixed.
- <br>
+ <br>
(WB, 2005/03/03)
</p>
sometimes declared to iterators equal if only they had the same
triangulation, refinement level, and index. Comparing iterators
into different containers is now completely disallowed.
- <br>
+ <br>
(WB, 2005/08/08)
</p>
vector-valued. It then makes up names for the individual components by
adding numbers to the name, just like the <code
class="class">DataOut</code> class already does for a long time.
- <br>
+ <br>
(WB, 2005/07/29)
</p>
New: The <code>DataOutStack</code> class can now also be
used to stack two-dimensional time or parameter dependent data into a 3d
output.
- <br>
+ <br>
(WB, 2005/07/25)
</p>
assumption that the coarse grid space is embedded into the fine
grid space. In particular, no interpolation points are
required.<p>
-
+
<p>Using this function, constraint matrices can be computed in
a general way.
- <br>
+ <br>
(GK, 2005/07/14)
</p>
<br>
(GK 2005/06/29)
</p>
-
+
<li> <p>
New: Added function <code
class="class">FiniteElementBase::component_to_system_index</code>
also for the case of non-standard oriented faces in 3d. This
avoids several awful pieces of code including questioning
face_orientation and using child_switch_tables.
- <br>
+ <br>
(RH, 2005/06/16)
</p>
a <code>FE_Q</code> object. Now this function can
also be called by the <code>FE_Q</code>
constructor which avoids code duplication.
- <br>
+ <br>
(RH, 2005/06/13)
</p>
<li> <p>
New: The method <code>create_mass_matrix</code>
- in class <code>MatrixCreator</code> can now handle
+ in class <code>MatrixCreator</code> can now handle
vector valued finite elements. A similar change was applied
to <code>create_right_hand_side</code> and
<code>create_boundary_right_hand_side</code>
now make the <code>project</code> function work also
for Raviart-Thomas elements and other vector valued FEs. This
is very useful, if some initial conditions have to be specified.
- <br>
+ <br>
(Oliver Kayser-Herold, 2005/06/03)
</p>
now shows the solution also on the deformed mesh. This
is an option and requires the mapping to be specified as
additional parameter to <code>build_patches</code>.
- <br>
+ <br>
(Oliver Kayser-Herold, 2005/05/31)
</p>
Gauss quadrature points on edges and in the interior for its
node values. The implementation is restricted to Cartesian mesh
cells right now, but works in 2D and 3D.
- <br>
+ <br>
(GK, 2005/05/30)
</p>
Improved: The <code>Mapping::transform_*</code>
functions accept <code>VectorSlice</code> instead
of <code>Vector</code>, thus allowing more flexibility.
- <br>
+ <br>
(GK, 2005/05/24)
</p>
<li> <p>
New: The <code>MatrixTools::apply_boundary_values</code>
function now also works for PETSc sequential and parallel matrices.
- <br>
+ <br>
(WB, 2005/05/05)
</p>
mechanism provided by the class <code
class="class">PathSearch</code>. Furthermore, a library of input
meshes has been started in <code>lib/meshes</code>.
- <br>
+ <br>
(GK, 2005/05/03)
</p>
overloaded and cell data was to be output; in that case, data from the
wrong cells was written out. This is now fixed. In contrast to this,
nodal data was never affected.
- <br>
+ <br>
(WB, 2005/04/20)
</p>
New: The new <code>GeometryInfo</code>::<code
class="member">line_to_cell_vertices</code> function maps line
vertex numbers to cell vertex numbers.
- <br>
+ <br>
(RH, 2005/03/11)
</p>
New: The new <code>GeometryInfo</code>::<code
class="member">face_to_cell_lines</code> function maps face
line numbers to cell line numbers.
- <br>
+ <br>
(RH, 2005/03/11)
</p>
New: The new <code>GeometryInfo</code>::<code
class="member">face_to_cell_vertices</code> function maps face
vertex numbers to cell vertex numbers.
- <br>
+ <br>
(RH, 2005/03/11)
</p>
class="class">KellyErrorEstimator</code> class that resulted in
assertions being thrown when run with multithreading
enabled. This is now fixed.
- <br>
+ <br>
(WB, 2005/03/10)
</p>
vertices. This leads to a change in the ordering of vertices in
output files generated by <code>GridOut</code>
for refined grids.
- <br>
+ <br>
(RH, 2005/03/09)
</p>
that happens when one coarsens an existing triangulation. In
that case, it would throw unjustified exceptions. This is now
fixed.
- <br>
+ <br>
(WB, 2005/03/01)
</p>
virtual abstract functions but no virtual destructor. This caused
warnings with some compilers, and is generally bad practice
anyway. This is now fixed.
- <br>
+ <br>
(WB, 2005/02/22)
</p>
automatic computation of embedding matrices under the sole
assumption that the coarse grid space is embedded into the fine
grid space. In particular, no interpolation points are required.
- <br>
+ <br>
(GK, 2005/02/08)
</p>
<li> <p>
Fixed: Several wrong assertions in the Raviart-Thomas finite element
class have been fixed, that were triggered when integrating face terms.
- <br>
- (Oliver Kayser-Herold, 2005/01/24;
+ <br>
+ (Oliver Kayser-Herold, 2005/01/24;
WB, 2005/01/31)
</p>
</ol>
namespace dealii;</code> declaration at the top of each file,
or qualify each deal.II class and function explicitly with
<code>dealii::</code> namespace access.
- <br>
+ <br>
(WB 2006/10/22)
</p>
Changed: Indices, such as vertex indices, are usually represented by
unsigned integers in deal.II. One place where we didn't do this was in
the <code>CellData</code> structure that can be used to describe cells
- when building a triangulation. This has now been rectified.
- <br>
+ when building a triangulation. This has now been rectified.
+ <br>
(WB 2006/09/06)
</p>
elements, the ordering of boundary DoFs returned by <code>DoFTools::map_dof_to_boundary_indices</code> has been
changed. Fortunately, this is a rarely used function so that the effect
should be limited to only a few programs.
- <br>
+ <br>
(WB 2006/04/26)
</p>
that you can't call these functions for face or edge iterators
any more, but only cells. Since this is what is usually
desired, this should not be too severe a restriction.
- <br>
+ <br>
(WB 2006/04/26)
</p>
Changed: The template argument of the <code>InterGridMap</code> class has been changed. Code
like e.g. <code>InterGridMap<DoFHandler,2></code> must be
replaced by <code>InterGridMap<DoFHandler<2> ></code>.
- <br>
+ <br>
(RH 2006/02/27)
</p>
argument is replaced by the type of DoFHandler they correspond
to. This may be <code>::DoFHandler</code> or
<code>hp::DoFHandler</code>.
- <br>
+ <br>
(RH 2006/02/10, WB 2005/12/20)
</p>
Removed: Support for gcc2.95 in particular, and all compilers that do
not support the <code>std::ostringstream</code> class in general,
has been removed.
- <br>
+ <br>
(WB 2006/02/03)
</p>
problems to application programs. On the other hand, it moves common
names like <code>Line</code> and <code>Quad</code> out of the global
namespace.
- <br>
+ <br>
(WB 2006/01/13)
</p>
action; if you really want the old behavior back, cast the
argument to the base class before the call, like so:
<code>DoFRenumbering::component_wise(static_cast<DoFHandler<dim>&>(mg_dof_handler))</code>)
- <br>
+ <br>
(WB 2005/12/15)
</p>
old vertex ordering. However, as this function might not be
supported forever users are advised to change their code to the
new numbering scheme and to use the <code>Triangulation::create_triangulation</code>
- function.
- <br>
+ function.
+ <br>
(RH 2005/11/02)
</p>
</ol>
</p>
<li> <p>New: The step-27 tutorial program has been added. It solves a Laplace
- equation with hp finite elements
+ equation with hp finite elements
and shows how to set finite element degrees, assemble matrices on
cells with varying polynomial degrees, and how to compute a simple
criterion for estimating the local Sobolev smoothness of a function.
been extended by new flags <tt>face_flip</tt> and <tt>face_rotation</tt>,
denoting a face which is rotated against the 'standard face' by 180 and 90
degrees, respectively. Furthermore, there is a new flag
- <tt>line_orientation</tt> with obvious meaning.
+ <tt>line_orientation</tt> with obvious meaning.
<br> These flags have to be respected during creation and refinement of a
triangulation, when projecting quadrature points to faces or when the dof
indices on a cell are extracted. Furthermore, asking for vertices and
lines of cells is in some cases a bit more complicated. The same applies,
- for example,
+ for example,
to the extraction of the information, which child of a neighbor is behind
a given subface. However, this information is supplied by various
functions in <code>GeometryInfo</code>. As a rule-of-thumb:
libraries that would contain illegal instructions. This is now fixed.
<br>
(WB 2006/08/09)
- </p>
+ </p>
<li> <p> Changed: The default number of subdivisions in the <code>build_patches</code> functions of the DataOut classes is
now part of the parameters read by <code>parse_parameters</code>. This default value is used
library instead of the contributed.
<br>
(GK 2006/03/23)
- </p>
+ </p>
<li> <p>New: The step-20 tutorial program was added. It shows how to use
Raviart-Thomas elements to solve a mixed formulation of the Laplace
the API manual (and the other way round) and leads to a
generally much nicer output. In particular, formulas embedded
in the documentation of the programs are now properly displayed.
- <br>
+ <br>
(WB 2006/01/16)
</p>
Improved: The link to the manual now points directly to
Doxygen's module list, and almost all classes are now grouped
into modules that capture related functionality.
- <br>
+ <br>
(WB 2006/01/16)
</p>
<code>PETSC_ARCH</code> environment variables but do not wish
deal.II to be configured for PETSc use, you
should specify <code>--with-petsc=no</code> as a flag during
- configuration.
- <br>
+ configuration.
+ <br>
(WB 2006/01/16)
</p>
impossible to figure out if you didn't know what was going
on. This mechanism has now been robustified and should yield
better error messages.
- <br>
+ <br>
(WB 2005/11/23)
</p>
New: The dynamic libraries, which is the name of shared libs
under Apples OSX, are now supported on these platforms and
enabled by default. They should safe a lot of harddisk space.
- <br>
+ <br>
(Oliver Kayser-Herold 2005/10/20)
</p>
<li> <p>
Fixed: The <tt>Tecplot</tt> library (<tt>tecio.a</tt>) was
detected but not added to <tt>LIBS</tt>. This is now fixed.
- <br>
+ <br>
(RH 2005/09/23)
</p>
<tt>NetCDF</tt> installation directory can also be specified
through the <tt>--with-netcdf=/path/to/netcdf</tt> configure
option.
- <br>
+ <br>
(RH 2005/09/23)
</p>
Fixed: The PETSc libraries have been relinked each time
<tt>make</tt> was called. Now, the PETSc libraries will only be
relinked if necessary.
- <br>
+ <br>
(RH 2005/09/15)
</p>
<li> <p> New: There is a new class:
<code>Functions::FEFieldFunction</code> which is a Function
- interface to a finite element solution.
- <br>
+ interface to a finite element solution.
+ <br>
(Luca Heltai 2007/08/29)
-
+
<li> <p> Improved: <code>FunctionDerivative</code> is now
derived from <code>AutoDerivativeFunction</code> and implements
gradients as well, giving you automatic second derivatives of a function.
<code>n</code> is a constant known at compile time. It allows to
calculate powers efficiently at compile time, most often things like a
number to the power <code>dim</code>.
- <br>
+ <br>
(WB 2007/06/23)
</p>
be of no concern to users. However, this restriction was
previously undocumented, even if implied. The documentation for
this has now been updated.
- <br>
+ <br>
(WB 2007/02/08)
</p>
<li> <p> New: The new <code>Functions::Monomial</code> class implements
monomials as a function object.
- <br>
+ <br>
(WB 2006/12/15)
</p>
<li> <p> Fixed: If no substring is found with a width smaller than the given
threshold, the <code>Utilities::break_text_into_lines</code> function now
returns the smallest possible substring (larger than the threshold).
- <br>
+ <br>
(Tobias Leicht 2006/12/13)
</p>
tensors. There is now a new version of this function, that contracts a
tensor of rank three with a second one of rank two over given indices
<tt>index1</tt> and <tt>index2</tt> of the first and second tensor,
- respectively, resulting in a tensor of rank three.
+ respectively, resulting in a tensor of rank three.
<br>
(Tobias Leicht 2006/09/27)
</p>
desirable, for example here:
<pre>
<code>
- void f()
+ void f()
{
Triangulation tria;
DoFHandler *dh = new DoFHandler(tria);
distance of a given point to the unit cell, and
<code>project_to_unit_cell</code> returning the projection of a point
onto the unit cell. Also, a new member <code>vertex_to_face</code> allow
- to determine to which faces of a cell a vertex belongs.
+ to determine to which faces of a cell a vertex belongs.
<br>
(Ralf B. Schulz 2006/05/10)
</p>
(GK 2006/03/08)
</p>
- <li> <p>Fixed: Writing a denormal "NaN" through <code>LogStream</code> objects such as
+ <li> <p>Fixed: Writing a denormal "NaN" through <code>LogStream</code> objects such as
<code>deallog</code> erroneously printed zero, rather
than "nan". This is now fixed.
<br>
<li> <p>
New: There are new functions
- <code>Utilities::match_at_string_start</code> and
+ <code>Utilities::match_at_string_start</code> and
<code>Utilities::get_integer_at_position</code>.
<br>
(WB 2005/12/19)
Fixed: The <tt>tecplot_binary</tt> <code>OutputFormat</code> has been added to <code>DataOutBase::get_output_format_names</code>. Now
an exception will be raised if <code>write_tecplot_binary</code> is invoked without
specifying the filename through the <code>DataOutBase::TecplotFlags</code> interface.
- <br>
+ <br>
(RH 2005/09/23)
</p>
<code>Utilities::System</code> namespace. In case
of a code not running in parallel, they simply return 1 and 0,
respectively.
- <br>
+ <br>
(RH 2005/09/20)
</p>
</ol>
(WB 2007/02/09)
</p>
- <li> <p>Improved: A simple <code>print_formatted</code>
+ <li> <p>Improved: A simple <code>print_formatted</code>
function has been added to <code>SparseMatrixEZ</code>.
<br>
(Moritz Allmaras 2007/02/08)
</p>
- <li> <p>Fixed: The <code>SparseDirectMA27</code>
+ <li> <p>Fixed: The <code>SparseDirectMA27</code>
class works for symmetric matrices and had a check to make sure the
matrix is indeed symmetric. However, the check compared entries for
equality, something one should not do for floating point numbers. It now
(WB 2006/12/21)
</p>
- <li> <p>Fixed: The <code>FullMatrix::invert</code>
+ <li> <p>Fixed: The <code>FullMatrix::invert</code>
function would return wrong results if one tried to invert a
matrix of dimension smaller than 5 in situ. This is now fixed.
<br>
now estimate the condition number of the linear system using <code>TridiagonalMatrix</code> and LAPACK.
<br>
(GK 2006/09/06)
- </p>
+ </p>
<li> <p>New: the class <code>TridiagonalMatrix</code>
has been introduced together with its basic methods, including an
</p>
<li> <p> Improved: PETSc 2.3.1 (<i>build 16</i>) is now supported by the
- linear preconditioner classes. The new PETSc functions
- <code>PCFactorSetXXX</code> are used.<br>
+ linear preconditioner classes. The new PETSc functions
+ <code>PCFactorSetXXX</code> are used.<br>
(Florian Prill, 2006/08/04)
</p>
the transpose of the matrix it points to.
<br>
(GK 2006/07/07)
- </p>
+ </p>
<li> <p>New: There is now a function <code>FullMatrix::trace</code> that does what its name suggests
it does.
<br>
(WB 2006/06/16)
</p>
-
+
<li> <p>Improved: <code>PointerMatrixAux</code> now has a
default constructor and a function for setting the <code>VectorMemory</code> object.
<br>
(GK 2006/06/14)
</p>
-
+
<li> <p>Fixed: <code>FullMatrix::print</code> would yield a link error
when used with <code>std::ofstream</code> since an explicit
instantiation was missing. The function has now been moved to
the header file and made inline so that it is always visible,
whatever the type of the output stream is.
- <br>
+ <br>
(WB 2006/06/07)
</p>
-
+
<li> <p> Improved: The <code>SparseDirectUMFPACK</code> solver can now also be
used with sparse matrices with elements of type float, as well
as with block matrices with types float and double.
- <br>
+ <br>
(WB 2006/04/25)
</p>
<li> <p> New: The function <code>BlockSparsityPattern::row_length</code> adds up
the row lengths of the individual blocks of a block matrix for
a given row.
- <br>
+ <br>
(WB 2006/04/25)
</p>
-
+
<li> <p> New: There is a new class <code>IdentityMatrix</code> that represents an
efficient version of the matrix with ones on the diagonal and
zeros everywhere else. The main usefulness of this matrix lies
assignment operators and copy constructors that allow to
conveniently fill a (sparse, full) matrix with the identity
matrix.
- <br>
+ <br>
(WB 2006/04/24)
</p>
-
+
<li> <p> New: There are now assignment operators from <code>BlockVector</code> to <code>Vector</code> and back.
<br>
(WB 2006/03/30)
the linear solver fails.
<br>
(Stephan Kramer, WB 2006/03/15)
- </p>
+ </p>
<li> <p>
Changed: There are new <code>FullMatrix</code>::<code>equ</code> functions which assign <tt>this</tt>
function for three vectors.
<br>
(RH 2006/02/21)
- </p>
+ </p>
<li> <p>
Fixed: The <code>SolverMinRes</code> class did not work
This is now fixed.
<br>
(Luca Heltai 2006/02/19)
- </p>
+ </p>
<li> <p>
Changed: There are now new <code>FullMatrix</code>::<code>add</code> functions which add two and three
scaled matrices.
<br>
(RH 2006/02/16)
- </p>
+ </p>
<li> <p>
Changed: The matrix classes <code>FullMatrix</code>, <code>SparseMatrix</code>, <code>SparseMatrixEZ</code> and <code>BlockMatrixBase</code> now have an <code>add</code> function analogous to the <code>add</code> of vector classes. The old <code>add_scaled</code> functions are now deprecated.
<br>
(RH 2006/02/16)
- </p>
+ </p>
- <li> <p>
+ <li> <p>
Improved: <code>BlockMatrixArray</code>::<code>enter_aux</code> allows using matrices without
adding vector multiplication by using <code>PointerMatrixAux</code>.
<br>
(GK 2006/02/02)
- </p>
+ </p>
<li> <p>
New: The class <code>PointerMatrixAux</code> was
Improved: The <code>FilteredMatrix</code> class
was able to filter only the <code>SparseMatrix</code>
matrix class.
- A new (faster) version is now in its place with the same interface
- that is able to perform such a filter on anything that provides
+ A new (faster) version is now in its place with the same interface
+ that is able to perform such a filter on anything that provides
the usual matrix vector operations.
- <br>
+ <br>
(Luca Heltai 2006/01/09)
</p>
specifically use a different memory allocation class, have been
changed to not specify anything at all, and thus fall back to
the default.
- <br>
+ <br>
(WB 2005/11/23)
</p>
(i.e. with a zero integer, rather than a floating point number on the
right), in that that opened up a second conversion sequence and the
compiler complained about an ambiguity. This is now fixed.
- <br>
+ <br>
(WB 2005/10/17)
</p>
also the <code>Tvmult</code> method. This feature
allows the straightforward use of <code>SparseILU</code>
as a multigrid smoother within the framework.
- <br>
+ <br>
(Oliver Kayser-Herold 2005/09/15)
</p>
</ol>
<li> <p> New: There is a new class:
<code>Functions::ParsedFunction</code> which is friendly
wrapper to the <code>FunctionParser</code> class.
- <br>
+ <br>
(Luca Heltai 2007/08/29)
-
+
<li> <p>Fixed: the function
- <code>DataOut::build_patches</code>
+ <code>DataOut::build_patches</code>
had a quadratic algorithm when generatic cell-data (as opposed
to DoF data). This algorithm became a bottleneck when
generating output on meshes with large number of cells. This is
</p>
<li> <p>New: the function
- <code>DoFTools::get_active_fe_indices</code>
+ <code>DoFTools::get_active_fe_indices</code>
extracts for each cell the active finite element index used on it.
<br>
(WB 2007/08/07)
</p>
<li> <p>Extended: the function
- <code>Triangulation::copy_triangulation</code>
+ <code>Triangulation::copy_triangulation</code>
copies all members of a triangulation except for the list of
<code>RefinementListener</code>s. In most cases this is exactly the
intended behavior. However, if a RefinementListener should be copied to
</p>
<li> <p>New: the function
- <code>DoFTools::make_sparsity_pattern</code>
+ <code>DoFTools::make_sparsity_pattern</code>
now takes an optional constraint matrix argument that can be used to
avoid the subsequent call to condense the sparsity pattern. In effect,
the sparsity pattern is generated condensed right away. For problems in
</p>
<li> <p>New: the new function
- <code>ConstraintMatrix::add_entries_local_to_global</code>
+ <code>ConstraintMatrix::add_entries_local_to_global</code>
can be used to add entries to a matrix where entries that will appear if
the current entry corresponds to a constrained degree of freedom are
also added.
</p>
<li> <p>Fixed: the function
- <code>GridTools::find_cells_adjacent_to_vertex</code>
+ <code>GridTools::find_cells_adjacent_to_vertex</code>
was not detecting properly the coarse cells adjacent to
refined cells.
<br>
<li> <p>Fixed: the two tools
<code>DoFTools::count_dofs_per_component</code> and
<code>DoFTools::count_dofs_per_block</code> where changing the
- size of the destination vector. Consistently with (most of) the
- rest of the library, now the vectors are expected to be the
- right size before calling these functions.
+ size of the destination vector. Consistently with (most of) the
+ rest of the library, now the vectors are expected to be the
+ right size before calling these functions.
<br>
(Luca Heltai 2007/05/15)
</p>
(WB 2007/02/20)
</p>
- <li> <p>New: Added function <code>GridGenerator</code>::<code >hyper_cube_with_cylindrical_hole</code> that produces
+ <li> <p>New: Added function <code>GridGenerator</code>::<code >hyper_cube_with_cylindrical_hole</code> that produces
a square with a circular hole in the middle in 2d, and extrudes it
- along the z-direction between 0 and L.
+ along the z-direction between 0 and L.
<br>
(Luca Heltai 2007/02/15)
</p>
be properly visualized. The boundary indicator was added to the
material flag of the faces (which is ignored when reading back the
mesh) in order for the Gmsh reader to be able to display the boundary
- indicator in a proper way.
+ indicator in a proper way.
<br>
(Luca Heltai 2007/02/15)
</p>
</p>
<li> <p>Fixed: <code>MatrixCreator</code>::<code >create_mass_matrix</code> and <code>
- create_laplace_matrix</code> computed wrong values for the right
- hand sides. This has been fixed.
+ create_laplace_matrix</code> computed wrong values for the right
+ hand sides. This has been fixed.
<br>
(Moritz Allmaras 2007/02/08)
</p>
<li> <p>
New: There are now two new functions <code>GridGenerator</code>::<code>subdivided_hyper_rectangle</code> that produces
a non-uniformly subdivided rectangle, ideally suited for graded
- meshes. One of these functions is able to create meshes with holes.
+ meshes. One of these functions is able to create meshes with holes.
<br>
(Yaqi Wang 2006/11/15, 2006/12/29)
</p>
- <li> <p> Fixed: Corrected <code>clone</code> method
- of <code>FE_DGQ</code> class for non-equidistant
+ <li> <p> Fixed: Corrected <code>clone</code> method
+ of <code>FE_DGQ</code> class for non-equidistant
support points.
<br>
(Florian Prill 2006/10/31)
New: There is a new functions <code>GridTools</code>::<code>create_union_triangulation</code>
that generates a triangulation that contains the respectively
finest cells of two input triangulations.
- <br>
+ <br>
(WB 2006/10/23)
</p>
of constraints are resolved at the time the
<code>ConstraintMatrix::close()</code> function is called. The only
thing that is not allowed are cycles in such constraints.
- <br>
+ <br>
(WB 2006/09/19)
</p>
<li> <p>
New: There are new functions <code>GridTools::minimal_cell_diameter</code> and <code>GridTools::maximal_cell_diameter</code>, with obvious
- functionality.
- <br>
+ functionality.
+ <br>
(WB 2006/09/06)
</p>
in particular did not allow to check for the number of matrices
in the array. The functions now take a reference to an array of
the correct length.
- <br>
+ <br>
(WB 2006/08/14)
</p>
<li> <p>
Extended: The <code>VectorTools::project</code> functions
are now also implemented for 1d.
- <br>
+ <br>
(WB 2006/08/08)
</p>
norm of this tensor is required later on, the new <code >DerivativeApproximation</code><code>::derivative_norm</code> function can be used. Note, that
for second order derivatives, this returns the largest eigenvalue
instead of the Frobenius norm, which is returned by the <code >Tensor<rank_,dim></code><code>::norm</code> function.
- <br>
+ <br>
(Tobias Leicht 2006/08/03)
</p>
symmetrized instead of an intermediate tensor. This should improve the
results slightly, but cause no problems otherwise, as this is completely
internal to the class.
- <br>
+ <br>
(Tobias Leicht 2006/08/03)
</p>
zero as the second argument (denoting the number of threads to use if
multithreading is enabled). This led to no output being created at
all. This is now fixed by throwing an exception in this case.
- <br>
+ <br>
(WB 2006/07/31)
</p>
<code>dofs_per_quad</code>, ..., depending on the explicitly
specified function template argument. This is often useful for
template trickery.
- <br>
+ <br>
(WB, 2006/07/28)
</p>
allowed (if all children are flagged for coarsening). This leeds to
triangulations with a slightly reduced number of cells. In some cases
older references will have to be updated.
- <br>
+ <br>
(Tobias Leicht 2006/06/22)
</p>
of a Dirac pulse on the base functions. This models a point source as
used in the calculations of Green's functions and can also be used to
extract the value of a finite element solution in a single point.
- <br>
+ <br>
(Ralf B. Schulz, 2006/05/15)
</p>
Changed: Functions <code>VectorTools</code>::<code>point_value</code> and <code>VectorTools</code>::<code>point_difference</code> using the old interface
without boundary mapping were replaced by wrapper functions
calling the new versions.
- <br>
+ <br>
(Ralf B. Schulz, 2006/05/15)
</p>
version as it automatically delivers also the local coordinate of the
point (so it can save some computation time in most cases as you don't have
to calculate that again).
- <br>
+ <br>
(Ralf B. Schulz, 2006/05/12)
</p>
<li> <p>
Improved: The functions <code>VectorTools</code>::<code>point_value</code> and <code>VectorTools</code>::<code>point_difference</code> now can also use arbitrary
mappings, using the new <code>GridTools</code>::<code>find_active_cell_around_point</code> algorithm.
- <br>
+ <br>
(Ralf B. Schulz, 2006/05/11)
</p>
determine the active cell in which a given point is located. For
points located on boundaries and edges, it is in most cases also able
to give the finest cell.
- <br>
+ <br>
(Ralf B. Schulz 2006/05/10)
</p>
moved to the <code>DoFCellAccessor</code> class that provides
access to cells (i.e. lines in 1d, quads in 2d, and hexes in
3d) for which this operation is actually useful.
- <br>
+ <br>
(WB 2006/05/01)
</p>
-
+
<li> <p> Fixed: second derivatives where not computed correctly in <code>FEFaceValuesBase</code>, i.e. when evaluating
second derivatives on faces of cells. This is now fixed. Using
second derivatives evaluated at quadrature points within a cell
<code>Triangulation::clear_user_flags_quad</code>, and
<code>Triangulation::clear_user_flags_hex</code> can be used to
selectively clear only some of the user flags as needed.
- <br>
+ <br>
(WB 2006/04/25)
</p>
New: The function <code>VectorTools::project</code> functions can now
also be used for vector arguments of type other than
<code>Vector<double></code>.
- <br>
+ <br>
(WB 2006/04/17)
</p>
to integrate the solution on one mesh against the shape functions on
another mesh, a situation that frequently happens in time-dependent but
also in nonlinear problems.
- <br>
+ <br>
(WB 2006/04/14)
</p>
New: The function <code>GridTools::have_same_coarse_mesh</code> figures
out whether two triangulations, two DoFHandlers, etc, are built on the
same coarse mesh.
- <br>
+ <br>
(WB 2006/04/14)
</p>
fairly frequent operation (called about 7 times per cell in
step-14, but much more often in more complicated programs), the
results of this function are now cached for faster access.
- <br>
+ <br>
(WB 2006/03/29)
</p>
element that also has non-primitive base elements (for example,
trying to interpolate boundary values for a the last component
of a RT x Q1 element). This is now fixed.
- <br>
+ <br>
(WB 2006/03/27)
</p>
New: There is a new function <code>FiniteElement</code>::<code>face_to_equivalent_cell_index</code> that can
convert a face index of a degree of freedom into a
corresponding cell index.
- <br>
+ <br>
(WB 2006/03/27)
</p>
<li> <p>
New: There are now functions <code>VectorTools</code>::<code>point_value</code> that evaluate the value of a
finite element function at a given point inside the domain.
- <br>
+ <br>
(WB 2006/03/06)
</p>
Improved: The <code>Triangulation</code>, <code>PersistentTriangulation</code>, <code>DoFHandler</code>, <code>hp::DoFHandler</code> and <code>MGDoFHandler</code> classes now all have a <code>dimension</code> variable which allows to ask
the template argument <tt>dim</tt> by
<code>SomeClass::dimension</code>.
- <br>
+ <br>
(RH 2006/02/27)
</p>
<li> <p>
Fixed: When used, the copy constructor of <code>MappingQ</code> would lead to memory
corruption. This is now fixed.
- <br>
+ <br>
(RH 2006/02/23)
</p>
</p>
<li> <p>
- New: There is now a <code>GridOut</code>::<code>write_gmsh</code> function to write in the gmsh
- format. It has the same features of the write_ucd one.
- <br>
+ New: There is now a <code>GridOut</code>::<code>write_gmsh</code> function to write in the gmsh
+ format. It has the same features of the write_ucd one.
+ <br>
(Luca Heltai 2006/02/17)
</p>
the standard <code>operator[]</code>
function. Similarly, the <code>FECollection</code>::<code>n_finite_elements</code>, <code>QCollection</code>::<code>n_quadratures</code> and <code>MappingCollection</code>::<code>n_mappings</code> functions are now renamed to
the standard <code>size()</code> function.
- <br>
+ <br>
(RH 2006/02/10)
</p>
<li> <p>
- Improved: <code>GridGenerator</code>::<code>half_hyper_ball</code> now is implemented also
- in three dimensions, and the boundary now is colored with 0 and one
+ Improved: <code>GridGenerator</code>::<code>half_hyper_ball</code> now is implemented also
+ in three dimensions, and the boundary now is colored with 0 and one
respectively on the curved and plane.
- <br>
+ <br>
(Luca Heltai 2006/02/07)
</p>
Fixed: The <code>MGDoFHandler</code>::<code>distribute_dofs</code> and <code>renumber</code> functions could not handle
coarsened grids (unused vertices, faces and cells). This is now
fixed.
- <br>
+ <br>
(RH 2006/01/23)
</p>
New: The <code>Mapping</code> class and derived classes now have
functions <code>clone()</code> that return a new object of the same
(derived) type.
- <br>
+ <br>
(WB 2006/1/8)
</p>
Improved: The <code>KellyErrorEstimator</code> class now also allows to
select the cells on which it is supposed to work by giving a material
ID, instead of only a subdomain ID.
- <br>
+ <br>
(WB 2006/1/4)
</p>
Improved: A new <code>TriaAccessor::ExcCellHasNoChildren</code>
exception will be raised if the <code>TriaObjectAccessor::child_index</code> function
is invoked for cells without children.
- <br>
+ <br>
(RH 2005/12/09)
</p>
thrown. In particular, this happened when the step-14 program was run on
dual processor machines and the library was compiled for
multithreading. This is now fixed.
- <br>
+ <br>
(WB 2005/10/20)
</p>
<li> <p>
New: There is now a new <code>Triangulation::get_boundary_indicators</code>
function.
- <br>
+ <br>
(RH 2005/09/30)
</p>
New: There is now a new <code>GridTools::delete_unused_vertices</code>
function. Previously a <tt>private</tt> function in <code>GridIn</code> it has now been moved to and made
<tt>public</tt> in <code>GridTools</code>.
- <br>
+ <br>
(RH 2005/09/28)
</p>
&filename)</code> function reads grids from NetCDF files. The
only data format currently supported is the <tt>TAU grid
format</tt>.
- <br>
+ <br>
(RH 2005/09/23)
</p>
function ran into an exception when called with
<tt>Default</tt> format and a filename which does not include
<tt>/</tt>. This is now fixed.
- <br>
+ <br>
(RH 2005/09/21)
</p>
</ol>
<br>
(WB, 2010/12/20)
</p></li>
-
+
<li><p>New: When computing errors using
VectorTools::integrate_difference in codimension one problems, if
you specified a norm that requires the computation of the gradients,
<br>
(Luca Heltai, 2010/12/05)
</p></li>
-
+
<li><p>New: The VectorTools::interpolate_boundary_values function can now
also be used to interpolate boundary values on meshes embedded in a higher
dimensional space.
<ol>
<li> New: VectorTools::get_position_vector now works with arbitrary
- FESystems, provided that the geometrical components are primitive,
- and that you provide a component mask to select what components of
+ FESystems, provided that the geometrical components are primitive,
+ and that you provide a component mask to select what components of
the finite element to use for the geometrical interpolation.
<br>
(Luca Heltai, 2015/07/25)
(David Wells, 2016/07/07)
</li>
- <li> Fixed: Fix various warnings introduced by GCC 6
+ <li> Fixed: Fix various warnings introduced by GCC 6
<br>
(David Wells, 2016/05/07)
</li>
</li>
<li>
- Fixed: The GridOut::write_mesh_per_processor_as_vtu() function now works for a mesh whose multilevel hierarchy is not distributed, as well as a mesh whose level_subdomain_ids do not necessarily match its subdomain_ids for every active cell.
+ Fixed: The GridOut::write_mesh_per_processor_as_vtu() function now works for a mesh whose multilevel hierarchy is not distributed, as well as a mesh whose level_subdomain_ids do not necessarily match its subdomain_ids for every active cell.
<br>
(Conrad Clevenger, 2017/02/09)
</li>
<li>
- Fixed: Patterns::Selection now ignores spaces at the beginning and the end
+ Fixed: Patterns::Selection now ignores spaces at the beginning and the end
of the sequence string.
<br>
(Rajat Arora, 2017/02/09)
<li>
Fixed: A bug in the rotation matrix used in step-18 has now been corrected.
This tutorial now uses the function
- Physics::Transformations::Rotations::rotation_matrix_3d to compute the 3d
+ Physics::Transformations::Rotations::rotation_matrix_3d to compute the 3d
rotation matrix.
<br>
(Jean-Paul Pelteret, Paul Kuberry, 2017/01/09)
</li>
<li>
- Improved: The MPI_Comm used in MatrixFree is not configurable via
+ Improved: The MPI_Comm used in MatrixFree is not configurable via
AdditionalData anymore but set to the MPI_Comm of the Triangulation.
<br>
(Daniel Arndt, 2017/01/08)
</li>
<li>
- New: Automatically recreate changes.h from files in
+ New: Automatically recreate changes.h from files in
subfolders of ./doc/news
<br>
(Daniel Arndt, 2016/11/18)
<li>
Fixed: Double contraction of two SymmetricTensor<..,VectorizedArray<T>> now
- works. Introduced internal::NumberType<T> with static member
+ works. Introduced internal::NumberType<T> with static member
internal::Numbertype::value to be called instead of the constructor Number()
in symmetric_tensor.h.
<br>
<ol>
<li>
- New: A new tutorial program, step-60, shows how to deal with non-matching discretizations using
- Distributed Lagrange Multipliers.
+ New: A new tutorial program, step-60, shows how to deal with non-matching discretizations using
+ Distributed Lagrange Multipliers.
<br>
(Luca Heltai, Giovanni Alzetta, 2018/05/07)
</li>
COMMAND ${CMAKE_COMMAND} -DOUTPUT_FILE=${CMAKE_CURRENT_BINARY_DIR}/changes.h
-P ${CMAKE_CURRENT_SOURCE_DIR}/changes/create_changes_h.cmake
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/changes
- DEPENDS
+ DEPENDS
${CMAKE_CURRENT_SOURCE_DIR}/changes/header
${CMAKE_CURRENT_SOURCE_DIR}/changes/header_incompatibilities
${CMAKE_CURRENT_SOURCE_DIR}/changes/header_major
<ul>
<li> The part "<code>/path/to/install/dir</code>" in the
command above needs to be replaced by the path to the directory which deal.II should be installed into. This can be a directory in your home directory (e.g., <code>~/bin/deal.II</code>) or a directory such as <code>/usr/local</code> if you have root privileges.
- Another option is to use something like <code>`pwd`/../installed/</code> (note the backticks).
+ Another option is to use something like <code>`pwd`/../installed/</code> (note the backticks).
</li>
<li> <code>make test</code> automatically runs in parallel and no <code>-jN</code> flag should be supplied to it.
while providing a multitude of linear solvers and preconditioners
and extension to other linear operators with minimal changes required in the code.
To enable Ginkgo, pass <code>-DGINKGO_DIR=/path/to/ginkgo</code> to CMake when
- configuring deal.II. For more detailed instructions, one can refer to
+ configuring deal.II. For more detailed instructions, one can refer to
<a href="external-libs/ginkgo.html" target="body">this page</a>.
</p>
</dd>
<code>-DSUNDIALS_DIR=/path/to/sundials</code> to the deal.II CMake call.
</p>
</dd>
-
+
<dt><a name="SymEngine"/>
<a href="https://symengine.github.io/" target="_top">SymEngine</a></dt>
<dd>
#
# Set up all executables:
#
- FILE(GLOB _steps
+ FILE(GLOB _steps
${CMAKE_CURRENT_SOURCE_DIR}/step-*/step-*.cc
${CMAKE_CURRENT_SOURCE_DIR}/step-*/step-*.cu)
FOREACH(_step ${_steps})
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
<h1>Results</h1>
-Running the program produces graphics of two grids (grid-1.svg and grid-2.svg).
+Running the program produces graphics of two grids (grid-1.svg and grid-2.svg).
You can open these with most every web browser -- in the simplest case,
just open the current directory in your file system explorer and click
on the file. If you like working on the command line, you call your
see any other way to find a bug, here's the one piece of
advice we will provide in this program: learn how to use a debugger as soon as
possible. It will be time well invested.
-(@dealiiVideoLectureSeeAlso{25}) The deal.II Frequently Asked
+(@dealiiVideoLectureSeeAlso{25}) The deal.II Frequently Asked
Questions (FAQ) page linked to from the top-level <a
href="http://www.dealii.org/">deal.II webpage</a> also provides a good number
of hints on debugging deal.II programs.
equation, while the second gives many examples of applications of these
concepts to a large number of other, more complicated equations. For
applications to individual types of equations, see also the publications by
-Becker @cite Bec95, @cite Bec98, Kanschat @cite Kan96, @cite FK97, Suttmeier
-@cite Sut96, @cite RS97, @cite RS98c, @cite RS99, Bangerth @cite BR99b,
+Becker @cite Bec95, @cite Bec98, Kanschat @cite Kan96, @cite FK97, Suttmeier
+@cite Sut96, @cite RS97, @cite RS98c, @cite RS99, Bangerth @cite BR99b,
@cite Ban00w, @cite BR01a, @cite Ban02, and Hartmann @cite Har02, @cite HH01,
-@cite HH01b. All of these works, from the original introduction by Becker and
-Rannacher to individual contributions to particular equations, have later been
-summarized in a book by Bangerth and Rannacher that covers all of these topics,
+@cite HH01b. All of these works, from the original introduction by Becker and
+Rannacher to individual contributions to particular equations, have later been
+summarized in a book by Bangerth and Rannacher that covers all of these topics,
see @cite BR03.
and so is $A$. We can thus use the CG method for solving the Newton steps.
(The fact that the matrix $A$ is symmetric and positive definite should not come
as a surprise. It results from taking the derivative of an operator that
-results from taking the derivative of an energy functional: the minimal
+results from taking the derivative of an energy functional: the minimal
surface equation simply minimizes some non-quadratic energy. Consequently,
the Newton matrix, as the matrix of second derivatives of a scalar energy,
must be symmetric since the derivative with regard to the $i$th and $j$th
degree of freedom should clearly commute. Likewise, if the energy functional
is convex, then the matrix of second derivatives must be positive definite,
-and the direct calculation above simply reaffirms this.)
+and the direct calculation above simply reaffirms this.)
It is worth noting, however, that the positive definiteness degenerates for
problems where $\nabla u$ becomes large. In other words, if we simply multiply
every processor in the MPI network only stores
a part of the matrix or vector. More specifically, each processor will
only store those rows of the matrix that correspond to a degree of
-freedom it "owns". For vectors, they either store only elements that
+freedom it "owns". For vectors, they either store only elements that
correspond to degrees of freedom the processor owns (this is what is
necessary for the right hand side), or also some additional elements
that make sure that every processor has access the solution components
-that live on the cells the processor owns (so-called
+that live on the cells the processor owns (so-called
@ref GlossLocallyActiveDof "locally active DoFs") or also on neighboring cells
(so-called @ref GlossLocallyRelevantDof "locally relevant DoFs").
Developing software to run in %parallel via MPI requires a bit of a change in
mindset because one typically has to split up all data structures so that
every processor only stores a piece of the entire problem. As a consequence,
-you can't typically access all components of a solution vector on each
+you can't typically access all components of a solution vector on each
processor -- each processor may simply not have enough memory to hold the
entire solution vector. Because data is split up or "distributed" across
processors, we call the programming model used by MPI "distributed memory
that multiple processors can all access all data within one memory
space, for example whenever multiple cores in a single machine work
on a common task). Some of the fundamentals of distributed memory
-computing are discussed in the
+computing are discussed in the
@ref distributed "Parallel computing with multiple processors using distributed memory"
documentation module, which is itself a sub-module of the
-@ref Parallel "Parallel computing" module.
+@ref Parallel "Parallel computing" module.
In general, to be truly able to scale to large numbers of processors, one
needs to split between the available processors <i>every</i> data structure
whose size scales with the size of the overall problem. (For a definition
-of what it means for a program to "scale", see
+of what it means for a program to "scale", see
@ref GlossParallelScaling "this glossary entry".) This includes, for
example, the triangulation, the matrix, and all global vectors (solution, right
hand side). If one doesn't split all of these objects, one of those will be
and not paying attention, then you're out of luck: you'll simply have to
wait until your requested over there will be worked on. In some cases,
bugs will lead the original receiver to never check your mail, and in that
-case you will wait forever -- this is called a <i>deadlock</i>.
+case you will wait forever -- this is called a <i>deadlock</i>.
(@dealiiVideoLectureSeeAlso{39,41,41.25,41.5})
In practice, one does not usually program at the level of sending and
The techniques this program then demonstrates are:
- How to use the PETSc wrapper classes; this will already be visible in the
declaration of the principal class of this program, <code>ElasticProblem</code>.
-- How to partition the mesh into subdomains; this happens in the
+- How to partition the mesh into subdomains; this happens in the
<code>ElasticProblem::setup_system()</code> function.
- How to parallelize operations for jobs running on an MPI network; here, this
is something one has to pay attention to in a number of places, most
information for the last step, though, but you get the idea. All this is after
release mode has been enabled by running <code>make release</code>, and
with the generation of graphical output switched off for the reasons stated in
-the program comments above.
+the program comments above.
(@dealiiVideoLectureSeeAlso{18})
The biggest 2d computations I did had roughly 7.1
million unknowns, and were done on 32 processes. It took about 40 minutes.
The program keeps a complete copy of the Triangulation and DoFHandler objects
on every processor. It also creates complete copies of the solution vector,
and it creates output on only one processor. All of this is obviously
-the bottleneck as far as parallelization is concerned.
+the bottleneck as far as parallelization is concerned.
Internally, within deal.II, parallelizing the data
structures used in hierarchic and unstructured triangulations is a hard
(\mathbf{f}, \varphi)_{\Omega(t_{n-1})}
-(\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}
\\
- &\qquad
+ &\qquad
+(\mathbf{b}(\mathbf{x},t_n)-\mathbf{b}(\mathbf{x},t_{n-1}), \varphi)_{\Gamma_N}
+(\sigma^{n-1} \mathbf{n}, \varphi)_{\Gamma_N}
\\
that reference all files for a given timestep. Sadly, the record files have a
different format between VisIt and Paraview, so we write out both formats.
-The code will generate the files <code>solution-TTTT.NNN.vtu</code>,
-where <code>TTTT</code> is the timestep number (starting from 1) and
+The code will generate the files <code>solution-TTTT.NNN.vtu</code>,
+where <code>TTTT</code> is the timestep number (starting from 1) and
<code>NNN</code> is the process rank (starting from
0). These files contain the locally owned cells for the timestep and
processor. The files <code>solution-TTTT.visit</code> is the visit record
you have to tell it about the @ref GlossMPICommunicator "MPI communicator"
that connects these processes upon construction of the triangulation).
Otherwise, the parallel::shared::Triangulation looks, for all practical
-purposes, like a regular Triangulation object.
+purposes, like a regular Triangulation object.
The convenience of using this class does not only result from being able
to avoid the manual call to GridTools::partition(). Rather, the DoFHandler
class now also knows that you want to use it in a parallel context, and
by default automatically enumerates degrees of freedom in such a way
that all DoFs owned by process zero come before all DoFs owned by process 1,
-etc. In other words, you can also avoid the call to
+etc. In other words, you can also avoid the call to
DoFRenumbering::subdomain_wise().
There are other benefits. For example, because the triangulation knows that
@endcode
This knowledge extends to the DoFHandler object built on such triangulations,
which can then identify which degrees of freedom are locally owned
-(see @ref GlossLocallyOwnedDofs) via calls such as
+(see @ref GlossLocallyOwnedDofs) via calls such as
DoFHandler::compute_n_locally_owned_dofs_per_processor() and
DoFTools::extract_locally_relevant_dofs(). Finally, the DataOut class
also knows how to deal with such triangulations and will simply skip
for a stable and accurate solution. The book and paper cited above give
indications as to how to do this, but there is also a large volume of
literature on this subject; a good start to get an overview of the topic can
-be found in the references of the paper by H.-Y. Duan and Q. Lin; @cite DL05.
+be found in the references of the paper by H.-Y. Duan and Q. Lin; @cite DL05.
<h5>Refinement during timesteps</h5>
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
and mesh. We will keep this step for the step-3 program and rather talk about
one practical aspect of a finite element program, namely that finite element
matrices are always very sparse: almost all entries in these
-matrices are zero.
+matrices are zero.
To be more precise, we say that a matrix is sparse
if the number of nonzero entries <i>per row</i> in the matrix is
problem of step-22 -- and in particular in 3d, the number of entries per row
may be several hundred. But the important point is that this number is
independent of the overall size of the problem: If you refine the mesh, the
-maximal number of unknowns per row remains the same.
+maximal number of unknowns per row remains the same.
Sparsity is one of the distinguishing feature of
the finite element method compared to, say, approximating the solution of a
${\cal O}(N)$ operations. It is clear that this wouldn't be possible if the
matrix were not sparse (because then the number of entries in the matrix would
have to be ${\cal O}(N^s)$ with some $s>1$, and doing a fixed number of
-matrix-vector products would take ${\cal O}(N^s)$ operations),
+matrix-vector products would take ${\cal O}(N^s)$ operations),
but it also requires very specialized solvers such as
multigrid methods to satisfy the requirement that the solution requires only a
fixed number of matrix-vector multiplications. We will frequently look at the
a DoFHandler, i.e., an object that represents all of the potentially nonzero
elements of a matrix one may build when discretizing a partial differential
equation on a mesh and its DoFHandler. This lack of structure in the sparsity
-pattern will be apparent from the pictures we show below.
+pattern will be apparent from the pictures we show below.
For most applications and algorithms, the exact way in which degrees of freedom
are numbered does not matter. For example, the Conjugate Gradient method we
On the other hand, if the fluid is confined in pores, then friction forces
exerted by the pore walls on the fluid become more and more important and
internal viscous friction becomes less and less important. Modeling this
-then first leads to the
+then first leads to the
<a href="https://en.wikipedia.org/wiki/Darcy%27s_law#Brinkman_form_of_Darcy's_law">Brinkman
model</a> if both effects are important, and in the limit of very small pores
to the <a href="https://en.wikipedia.org/wiki/Darcy%27s_law">Darcy equations</a>.
-The latter is just a different name for the Poisson or Laplace equation,
+The latter is just a different name for the Poisson or Laplace equation,
connotating it with the area to which one wants to apply it: slow flow
in a porous medium. In essence it says that the velocity is proportional
to the pressure gradient that drives the fluid through the porous medium.
pressure test function $q$), and it is indefinite. At least it is
symmetric. In other words: the Conjugate Gradient method is not going to
work since it is only applicable to problems in which the matrix is
-symmetric and positive definite.
+symmetric and positive definite.
We would have to resort to other iterative solvers instead, such as
MinRes, SymmLQ, or GMRES, that can deal with indefinite systems. However, then
the next problem immediately surfaces: due to the zero block, there are zeros
B^TM^{-1}B P &=& B^TM^{-1} F - G, \\
MU &=& F - BP.
@f}
-Here, the matrix $S=B^TM^{-1}B$ (called the
-<a href="https://en.wikipedia.org/wiki/Schur_complement">Schur complement</a>
+Here, the matrix $S=B^TM^{-1}B$ (called the
+<a href="https://en.wikipedia.org/wiki/Schur_complement">Schur complement</a>
of $A$)
is obviously symmetric and, owing to the positive definiteness of $M$ and the
fact that $B$ has full column rank, $S$ is also positive
method to it. However, computing $S$ is expensive because it requires us
to compute the inverse of the (possibly large) matrix $M$; and $S$ is in fact
also a full matrix because even though $M$ is sparse, its inverse $M^{-1}$
-will generally be a dense matrix.
+will generally be a dense matrix.
On the other hand, the CG algorithm doesn't require
us to actually have a representation of $S$: It is sufficient to form
matrix-vector products with it. We can do so in steps, using the fact that
matrix products are associative (i.e., we can set parentheses in such a
-way that the product is more convenient to compute):
+way that the product is more convenient to compute):
To compute $Sv=(B^TM^{-1}B)v=B^T(M^{-1}(Bv))$, we
<ol>
<li> compute $w = B v$;
matrix-vector products.
In the following, we will then have to come up with ways to represent the
-matrix $S$ so that it can be used in a Conjugate Gradient solver,
+matrix $S$ so that it can be used in a Conjugate Gradient solver,
as well as to define ways in which we can precondition the solution
of the linear system involving $S$, and deal with solving linear systems
-with the matrix $M$ (the second step above).
+with the matrix $M$ (the second step above).
@note The key point in this consideration is to recognize that to implement
an iterative solver such as CG or GMRES, we never actually need the actual
The second part of the equations is the description of the
dynamics of the saturation, i.e., how the relative concentration of the
-two fluids changes with time. The saturation equation for the displacing
+two fluids changes with time. The saturation equation for the displacing
fluid (water) is given by the following conservation law:
@f{eqnarray*}
S_{t} + \nabla \cdot (F(S) \mathbf{u}) = q_{w},
which can be rewritten by using the product rule of the divergence operator
in the previous equation:
@f{eqnarray*}
- S_{t} + F(S) \left[\nabla \cdot \mathbf{u}\right]
+ S_{t} + F(S) \left[\nabla \cdot \mathbf{u}\right]
+ \mathbf{u} \cdot \left[ \nabla F(S)\right]
= S_{t} + F(S) q + \mathbf{u} \cdot \nabla F(S) = q_{w}.
@f}
refined meshes and continuous pressure finite elements so that there
are hanging nodes. Another cause for such constraints are Dirichlet
boundary conditions on the pressure. Then the AffineConstraints
-class, upon copying the local contributions to the matrix into the
+class, upon copying the local contributions to the matrix into the
global linear system will zero out rows and columns corresponding
to constrained degrees of freedom and put a positive entry on
the diagonal. (You can think of this entry as being one for
Apart from this, many other algorithms have been tested and improved during
the creation of this program. For example, in building the sparsity pattern,
-we originally used a (now no longer existing) BlockCompressedSparsityPattern
+we originally used a (now no longer existing) BlockCompressedSparsityPattern
object that added one element at a time; however, its data structures were poorly
-adapted for the large numbers of nonzero entries per row created by our
-discretization in 3d, leading to a quadratic behavior. Replacing the internal
+adapted for the large numbers of nonzero entries per row created by our
+discretization in 3d, leading to a quadratic behavior. Replacing the internal
algorithms in deal.II to set many elements at a time, and using a
BlockCompressedSimpleSparsityPattern (which has, as of early 2015, been in turn
replaced by BlockDynamicSparsityPattern) as a better adapted data structure,
\end{array}\right)=\\
\left(\begin{array}{l}
G_1 \\
- G_2 -(\theta F^{n}+(1-\theta)F ^{n-1})c_{0}^{2}k \\
+ G_2 -(\theta F^{n}+(1-\theta)F ^{n-1})c_{0}^{2}k \\
\end{array}\right)
@f]
G_2 \\
\end{array} \right)=\\
\left(\begin{array}{l}
- M\bar{p}^{n-1}+k(1-\theta)Mv^{n-1}\\
+ M\bar{p}^{n-1}+k(1-\theta)Mv^{n-1}\\
(-c_{0}^{2}k (1-\theta)A+c_0 B)\bar{p}^{n-1} +Mv^{n-1}
\end{array}\right)
@f]
By simple transformations, one then obtains two equations for
the pressure potential and its derivative, just as in the previous tutorial program:
@f{eqnarray*}
-(M+(k\,\theta\,c_{0})^{2}A+c_0k\theta B)\bar{p}^{n} & = &
+(M+(k\,\theta\,c_{0})^{2}A+c_0k\theta B)\bar{p}^{n} & = &
G_{1}+(k\, \theta)G_{2}-(c_0k)^2\theta (\theta F^{n}+(1-\theta)F^{n-1}) \\
-Mv^n & = & -(c_0^2\,k\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 -
-c_0^2k(\theta F^{n}+(1-\theta)F^{n-1})
+Mv^n & = & -(c_0^2\,k\, \theta\, A+c_0B)\bar{p}^{n}+ G_2 -
+c_0^2k(\theta F^{n}+(1-\theta)F^{n-1})
@f}
&=
f(\mathbf x, t),
\qquad\qquad &&
- \forall \mathbf x \in \Omega, t\in (0,T),
+ \forall \mathbf x \in \Omega, t\in (0,T),
\\
u(\mathbf x, 0) &= u_0(\mathbf x) &&
\forall \mathbf x \in \Omega, \\
(1-\theta)f(\mathbf x, t_{n-1})
+
\theta f(\mathbf x, t_n)
- \right].
+ \right].
@f}
Here, $k_n=t_n-t_{n-1}$ is the time step size. The theta-scheme generalizes
the explicit Euler ($\theta=0$), implicit Euler ($\theta=1$) and
(1-\theta)f(\mathbf x, t_{n-1})
+
\theta f(\mathbf x, t_n)
- \right].
+ \right].
@f}
We can here consider $u^{n-1}$ as data since it has presumably been computed
- before. Now, let us replace
+ before. Now, let us replace
@f{align*}
- u^n(\mathbf x)\approx u_h^n(\mathbf x)
+ u^n(\mathbf x)\approx u_h^n(\mathbf x)
=
\sum_j U^n \varphi_j(\mathbf x),
@f}
@f}
Here,
@f{align*}
- \chi_1(\mathbf x) &=
+ \chi_1(\mathbf x) &=
\left\{
\begin{array}{ll}
1
\end{array}
\right.
\\
- \chi_2(\mathbf x) &=
+ \chi_2(\mathbf x) &=
\left\{
\begin{array}{ll}
1
If you look at the meshes in the movie above, it is clear that they are not
particularly well suited to the task at hand. In fact, they look rather
-random.
+random.
There are two factors at play. First, there are some islands where cells
have been refined but that are surrounded by non-refined cells (and there
but they also don't help because so many of their additional degrees of
freedom are in fact constrained by hanging node constraints. That said,
this is easy to fix: the Triangulation class takes an argument to its
-constructor indicating a level of "mesh smoothing". Passing one of many
-possible flags, this instructs the triangulation to refine some additional
-cells, or not to refine some cells, so that the resulting mesh does not have
+constructor indicating a level of "mesh smoothing". Passing one of many
+possible flags, this instructs the triangulation to refine some additional
+cells, or not to refine some cells, so that the resulting mesh does not have
these artifacts.
-The second problem is more severe: the mesh appears to lag the solution.
+The second problem is more severe: the mesh appears to lag the solution.
The underlying reason is that we only adapt the mesh once every fifth
time step, and only allow for a single refinement in these cases. Whenever a
source switches on, the solution had been very smooth in this area before and
SparseMatrixIterators::Iterator<double,true> row_it (&mass_matrix, 0);
-for(unsigned int m = 0; m<num_rows; ++m)
+for(unsigned int m = 0; m<num_rows; ++m)
{
// check the diagonal element
row_it = mass_matrix.begin(m);
mass_matrix_min_diag = std::min(row_it->value(), mass_matrix_min_diag);
++row_it;
-
+
// check the off-diagonal elements
for(; row_it != mass_matrix.end(m); ++row_it)
mass_matrix_max_offdiag = std::max(row_it->value(), mass_matrix_max_offdiag);
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
GridGenerator::hyper_cube(), as before:
@code
for (auto &face : triangulation.active_face_iterators())
- if (std::fabs(face->center()(1) - (-1.0)) < 1e-12 ||
+ if (std::fabs(face->center()(1) - (-1.0)) < 1e-12 ||
std::fabs(face->center()(1) - (1.0)) < 1e-12)
face->set_boundary_id(1);
@endcode
%HDF5 is a commonly used format that can be read by many scripting
languages (e.g. R or Python). It is not difficult to get deal.II to
produce some %HDF5 files that can then be used in external scripts to
-postprocess some of the data generated by this program. Here are some
+postprocess some of the data generated by this program. Here are some
ideas on what is possible.
<h4> Changing the output to .h5</h4>
-To fully make use of the automation we first need to introduce a private variable for the number of
-global refinement steps <code>unsigned int n_refinement_steps </code>, which will be used for the output filename.
+To fully make use of the automation we first need to introduce a private variable for the number of
+global refinement steps <code>unsigned int n_refinement_steps </code>, which will be used for the output filename.
In <code>make_grid()</code> we then replace <code>triangulation.refine_global(5);</code> with
@code
n_refinement_steps = 5;
The deal.II library has two different %HDF5 bindings, one in the HDF5
namespace (for interfacing to general-purpose data files)
and another one in DataOut (specifically for writing files for the
-visualization of solutions).
-Although the HDF5 deal.II binding supports both serial and MPI, the %HDF5 DataOut binding
+visualization of solutions).
+Although the HDF5 deal.II binding supports both serial and MPI, the %HDF5 DataOut binding
only supports parallel output.
For this reason we need to initialize an MPI
communicator with only one processor. This is done by adding the following code.
more general file format, it can also easily be processed in scripting
languages for other purposes.
-
+
<h4> Adding the point value and the mean (see extension above) into the .h5 file</h4>
After outputing the solution, the file can be opened again to include
information on the possible output options.)
To make this happen, we first include the necessary header into our file:
-@code
+@code
#include <deal.II/base/hdf5.h>
@endcode
-Adding the following lines to the end
+Adding the following lines to the end
of our output routine adds the information about the value of the
solution at a particular point, as well as the mean value of the
solution, to our %HDF5 file:
<a href="https://en.wikipedia.org/wiki/R_(programming_language)">R
programming language</a>, a widely used language in statistical data
analysis. (Similar things can also be done in Python, for example.)
-If you are unfamiliar with R and ggplot2 you could check out the data carpentry course on R
+If you are unfamiliar with R and ggplot2 you could check out the data carpentry course on R
<a href="https://datacarpentry.org/R-ecology-lesson/index.html">here</a>.
-Furthermore, since most search engines struggle with searches of the form "R + topic",
+Furthermore, since most search engines struggle with searches of the form "R + topic",
we recommend using the specializes service <a
href="http://rseek.org">RSeek </a> instead.
@endcode
This gives the following output
@code{.unparsed}
-HDF5 FILE
+HDF5 FILE
name /
-filename
+filename
name otype dclass dim
0 cells H5I_DATASET INTEGER x 1024
-1 mean_value H5I_DATASET FLOAT 1
+1 mean_value H5I_DATASET FLOAT 1
2 nodes H5I_DATASET FLOAT x 1089
-3 point_value H5I_DATASET FLOAT 1
+3 point_value H5I_DATASET FLOAT 1
4 solution H5I_DATASET FLOAT x 1089
@endcode
The datasets can be accessed by <code>h5f\$name</code>. The function
that is used to store our cells.
We can see the following three matrices, as well as the two
additional data points we added.
-<ul>
-<li> <code>cells</code>: a 4x1024 matrix that stores the (C++) vertex indices for each cell
+<ul>
+<li> <code>cells</code>: a 4x1024 matrix that stores the (C++) vertex indices for each cell
<li> <code>nodes</code>: a 2x1089 matrix storing the position values (x,y) for our cell vertices
<li> <code>solution</code>: a 1x1089 matrix storing the values of our solution at each vertex
</ul>
At first the data needs to be manipulated and added to a <code>data.frame</code>.
After that, a <code>ggplot</code> object is constructed and manipulated by adding plot elements to it.
-<code>nodes</code> and <code>cells</code> contain all the information we need to plot our grid.
+<code>nodes</code> and <code>cells</code> contain all the information we need to plot our grid.
The following code wraps all the data into one dataframe for plotting our grid:
@code{.r}
# Counting in R starts at 1 instead of 0, so we need to increment all
@code{.r}
pdf (paste("grid_",refinement,".pdf",sep=""),width = 5,height = 5) # Open new PDF file
plt <- ggplot(meshdata,aes(x=x,y=y,group=id)) # Construction of our plot
- # object, at first only data
+ # object, at first only data
plt <- plt + geom_polygon(fill="white",colour="black") # Actual plotting of the grid as polygons
plt <- plt + ggtitle(paste("grid at refinement level #",refinement))
more interesting.
To make a 2D pseudocolor plot of our solution we will use <code>geom_raster</code>.
This function needs a structured grid, i.e. uniform in x and y directions.
-Luckily our data at this point is structured in the right way.
+Luckily our data at this point is structured in the right way.
The following code plots a pseudocolor representation of our surface into a new PDF:
@code{.r}
pdf (paste("pseudocolor_",refinement,".pdf",sep=""),width = 5,height = 4.2) # Open new PDF file
</td>
</tr>
</table>
-<li>
+<li>
For plotting the converge curves we need to re-run the C++ code multiple times with different values for <code>n_refinement_steps</code>
-starting from 1.
+starting from 1.
Since every file only contains a single data point we need to loop over them and concatenate the results into a single vector.
@code{.r}
n_ref <- 8 # Maximum refinement level for which results are existing
<td><img src="https://www.dealii.org/images/steps/developer/step-3.extensions.convergence_mean.png" alt=""></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-3.extensions.convergence_point.png" alt=""></td>
</tr>
- </table>
+ </table>
</ul>
<h4>Cache the explicit part of residual</h4>
The residual calculated in ConservationLaw::assemble_cell_term function
-reads
+reads
$R_i = \left(\frac{\mathbf{w}^{k}_{n+1} - \mathbf{w}_n}{\delta t}
, \mathbf{z}_i \right)_K +
\theta \mathbf{B}({\mathbf{w}^{k}_{n+1}})(\mathbf{z}_i)_K +
fact that the solution of a pure Neumann problem is known up to an
arbitrary constant $c$, which means that, if we set the Neumann data
to be zero, then any constant $\phi = \phi_\infty$ will be a solution.
-Inserting the constant solution and the Neumann boundary condition in the
+Inserting the constant solution and the Neumann boundary condition in the
boundary integral equation, we have
@f{align*}
\alpha\left(\mathbf{x}\right)\phi\left(\mathbf{x}\right)
+\int_{\Gamma}\frac{ \partial G(\mathbf{y}-\mathbf{x}) }{\partial \mathbf{n}_y} \, ds_y
\right]
@f}
-The integral on $\Gamma_\infty$ is unity, see above, so division by the constant $\phi_\infty$ gives us the explicit
+The integral on $\Gamma_\infty$ is unity, see above, so division by the constant $\phi_\infty$ gives us the explicit
expression above for $\alpha(\mathbf{x})$.
While this example program is really only focused on the solution of the
<i>This program was contributed by Toby D. Young and Wolfgang
Bangerth. </i>
-<a name="Preamble"></a>
+<a name="Preamble"></a>
<h1>Preamble</h1>
The problem we want to solve in this example is an eigenspectrum
$\Psi_\ell,\varepsilon_\ell$, $\ell\in {\cal I}$, where $\cal I$ can
be a finite or infinite set (and in the latter case it may be a discrete or
sometimes at least in part a continuous set). In either case, let us note that
-there is
+there is
no longer just a single solution, but a set of solutions (the various
eigenfunctions and corresponding eigenvalues) that we want to
compute. The problem of numerically finding all eigenvalues
(nondimensionalized) form
@f{align*} [-\Delta + V(\mathbf x)]
\Psi(\mathbf x) &= \varepsilon \Psi(\mathbf x) \qquad &&\text{in}\
-\Omega\quad, \\ \Psi &= 0 &&\text{on}\ \partial\Omega\quad.
-@f}
+\Omega\quad, \\ \Psi &= 0 &&\text{on}\ \partial\Omega\quad.
+@f}
As a consequence, this particle can only exist in a certain number of
eigenstates that correspond to the energy eigenvalues
$\varepsilon_\ell$ admitted as solutions of this equation. The
finite element method gains much of its power and simplicity from
the fact that we just do the same thing on every cell, without
having to think too much about where a cell is, whether it bounds
-on a less refined cell and consequently has a hanging node, or is
+on a less refined cell and consequently has a hanging node, or is
adjacent to the boundary. All such checks would make the assembly
of finite element linear systems unbearably difficult to write and
-even more so to read.
+even more so to read.
Consequently, of course, when you distribute degrees of freedom with
your DoFHandler object, you don't care whether some of the degrees
of freedom you enumerate are at a Dirichlet boundary. They all get
numbers. We just have to take care of these degrees of freedom at a
later time when we apply boundary values. There are two basic ways
-of doing this (either using MatrixTools::apply_boundary_values()
+of doing this (either using MatrixTools::apply_boundary_values()
<i>after</i> assembling the linear system, or using
AffineConstraints::distribute_local_to_global() <i>during</i> assembly;
see the @ref constraints "constraints module" for more information),
but both result in the same: a linear system that has a total
number of rows equal to the number of <i>all</i> degrees of freedom,
-including those that lie on the boundary. However, degrees of
+including those that lie on the boundary. However, degrees of
freedom that are constrained by Dirichlet conditions are separated
from the rest of the linear system by zeroing out the corresponding
row and column, putting a single positive entry on the diagonal,
F_i \\ F_b
\end{pmatrix}.
@f}
-Here, subscripts $i$ and $b$ correspond to interior and boundary
+Here, subscripts $i$ and $b$ correspond to interior and boundary
degrees of freedom, respectively. The interior degrees of freedom
satisfy the linear system $A_i U_i=F_i$ which yields the correct
solution in the interior, and boundary values are determined by
from the process of eliminating boundary degrees of freedom, and
$F_b$ is chosen in such a way that $U_{b,j}=D_{b,jj}^{-1} F_{b,j}$
has the correct boundary values for every boundary degree of freedom
-$j$. (For the curious, the entries of the
+$j$. (For the curious, the entries of the
matrix $D_b$ result from adding modified local contributions to the
-global matrix where for the local matrices the diagonal elements, if non-zero,
-are set to their absolute value; otherwise, they are set to the average of
+global matrix where for the local matrices the diagonal elements, if non-zero,
+are set to their absolute value; otherwise, they are set to the average of
absolute values of the diagonal. This process guarantees that the entries
of $D_b$ are positive and of a size comparable to the rest of the diagonal
entries, ensuring that the resulting matrix does not incur unreasonable
we compute in this program happens to be in this interval, and in that
case we would not know immediately whether it is a spurious or a true
eigenvalue. In that case, one could simply scale the diagonal elements of
-either matrix after computing the two matrices,
+either matrix after computing the two matrices,
thus shifting them away from the frequency of interest in the eigen-spectrum.
This can be done by using the following code, making sure that all spurious
eigenvalues are exactly equal to $1.234\cdot 10^5$:
stiffness_matrix.set(i, i, 1.234e5);
mass_matrix.set(i, i, 1);
}
-@endcode
+@endcode
However, this strategy is not pursued here as the spurious eigenvalues
-we get from our program as-is happen to be greater than the lowest
+we get from our program as-is happen to be greater than the lowest
five that we will calculate and are interested in.
<h3>Implementation details</h3>
-The program below is essentially just a slightly modified version of
+The program below is essentially just a slightly modified version of
step-4. The things that are different are the following:
-<ul>
+<ul>
<li>The main class (named <code>EigenvalueProblem</code>) now no
longer has a single solution vector, but a whole set of vectors for
</ul>
-The rest of the program follows in a pretty straightforward way from
+The rest of the program follows in a pretty straightforward way from
step-4.
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
The program produces two files: <code>solution-2d.vtk</code> and
<code>solution-3d.vtk</code>, which can be viewed using the programs
Visit or Paraview (in case you do not have these programs, you can easily
-change the
+change the
output format in the program to something which you can view more
easily). Visualizing solutions is a bit of an art, but it can also be fun, so
you should play around with your favorite visualization tool to get familiar
axis, puts a big fat label "X Axis" on the $x$ axis and similar for the other
axes, shows the file name from which the data was taken in the top left and
the name of the user doing so and the time and date on the bottom right. None
-of this is important
+of this is important
here: the axes are equally easy to make out because the tripod at the bottom
left is still visible, and we know from the program that the domain is
$[-1,1]^3$, so there is no need for tick marks. As a consequence, I have
<li> The program allows you to select from two different coarse meshes
through the parameter file. These are either a cube $[0,1]^3$ or
a half sphere with the open side facing the positive $z$ direction.
-
+
<li>In either case, we will assume the convention that the part of the
boundary that may be in contact with the obstacle has boundary
indicator one. For both kinds of meshes, we assume that this is a free
Cycle 0:
Number of active cells: 512
Number of degrees of freedom: 2187
-
+
Newton iteration 1
Updating active set...
Size of active set: 1
- Assembling system...
- Solving system...
+ Assembling system...
+ Solving system...
Error: 173.076 -> 1.64265e-06 in 7 Bicgstab iterations.
Accepting Newton solution with residual: 1.64265e-06
-
+
Newton iteration 2
Updating active set...
Size of active set: 1
- Assembling system...
- Solving system...
+ Assembling system...
+ Solving system...
Error: 57.3622 -> 3.23721e-07 in 8 Bicgstab iterations.
Accepting Newton solution with residual: 24.9028
Active set did not change!
-
+
Newton iteration 3
Updating active set...
Size of active set: 1
- Assembling system...
- Solving system...
+ Assembling system...
+ Solving system...
Error: 24.9028 -> 9.94326e-08 in 7 Bicgstab iterations.
Residual of the non-contact part of the system: 1.63333
with a damping parameter alpha = 1
Newton iteration 6
Updating active set...
Size of active set: 1
- Assembling system...
- Solving system...
+ Assembling system...
+ Solving system...
Error: 1.43188e-07 -> 3.56218e-16 in 8 Bicgstab iterations.
Residual of the non-contact part of the system: 4.298e-14
with a damping parameter alpha = 1
Cycle 3:
Number of active cells: 14652
Number of degrees of freedom: 52497
-
+
Newton iteration 1
Updating active set...
Size of active set: 145
- Assembling system...
- Solving system...
+ Assembling system...
+ Solving system...
Error: 296.309 -> 2.72484e-06 in 10 Bicgstab iterations.
Accepting Newton solution with residual: 2.72484e-06
-
+
...
Newton iteration 10
Updating active set...
Size of active set: 145
- Assembling system...
- Solving system...
+ Assembling system...
+ Solving system...
Error: 2.71541e-07 -> 1.5428e-15 in 27 Bicgstab iterations.
Residual of the non-contact part of the system: 1.89261e-13
with a damping parameter alpha = 1
The picture shows the adaptive refinement and as well how much a cell is
plastified during the contact with the ball. Remember that we consider the
norm of the deviator part of the stress in each quadrature point to
-see if there is elastic or plastic behavior.
+see if there is elastic or plastic behavior.
The blue
color means that this cell contains only elastic quadrature points in
contrast to the red cells in which all quadrature points are plastified.
-In the middle of the top surface -
+In the middle of the top surface -
where the mesh is finest - a very close look shows the dimple caused by the
obstacle. This is the result of the <code>move_mesh()</code>
function. However, because the indentation of the obstacles we consider here
<i>This program was contributed by Jean-Paul Pelteret and Andrew McBride.
<br>
This material is based upon work supported by the German Science Foundation (Deutsche
-Forschungsgemeinschaft, DFG), grant STE 544/39-1, and the National Research Foundation of South Africa.
+Forschungsgemeinschaft, DFG), grant STE 544/39-1, and the National Research Foundation of South Africa.
</i>
@dealiiTutorialDOI{10.5281/zenodo.439772,https://zenodo.org/badge/DOI/10.5281/zenodo.439772.svg}
<a name="Intro"></a>
<h1>Introduction</h1>
-The subject of this tutorial is nonlinear solid mechanics.
-Classical single-field approaches (see e.g. step-18) can not correctly describe the response of quasi-incompressible materials.
+The subject of this tutorial is nonlinear solid mechanics.
+Classical single-field approaches (see e.g. step-18) can not correctly describe the response of quasi-incompressible materials.
The response is overly stiff; a phenomenon known as locking.
-Locking problems can be circumvented using a variety of alternative strategies.
+Locking problems can be circumvented using a variety of alternative strategies.
One such strategy is the three-field formulation.
It is used here to model the three-dimensional, fully-nonlinear (geometrical and material) response of an isotropic continuum body.
The material response is approximated as hyperelastic.
We begin with a crash-course in nonlinear kinematics.
For the sake of simplicity, we restrict our attention to the quasi-static problem.
Thereafter, various key stress measures are introduced and the constitutive model described.
-We then describe the three-field formulation in detail prior to explaining the structure of the class used to manage the material.
+We then describe the three-field formulation in detail prior to explaining the structure of the class used to manage the material.
The setup of the example problem is then presented.
@note This tutorial has been developed (and is described in the introduction) for the problem of elasticity in three dimensions.
While the space dimension could be changed in the main() routine, care needs to be taken.
Two-dimensional elasticity problems, in general, exist only as idealizations of three-dimensional ones.
- That is, they are either plane strain or plane stress.
- The assumptions that follow either of these choices needs to be consistently imposed.
+ That is, they are either plane strain or plane stress.
+ The assumptions that follow either of these choices needs to be consistently imposed.
For more information see the note in step-8.
<h3>List of references</h3>
The three-field formulation implemented here was pioneered by Simo et al. (1985) and is known as the mixed Jacobian-pressure formulation.
Important related contributions include those by Simo and Taylor (1991), and Miehe (1994).
-The notation adopted here draws heavily on the excellent overview of the theoretical aspects of nonlinear solid mechanics by Holzapfel (2001).
+The notation adopted here draws heavily on the excellent overview of the theoretical aspects of nonlinear solid mechanics by Holzapfel (2001).
A nice overview of issues pertaining to incompressible elasticity (at small strains) is given in Hughes (2000).
<ol>
pressure in solid mechanics as
$p = - 1/3 \textrm{tr} \boldsymbol{\sigma} = - 1/3 J^{-1} \textrm{tr} \boldsymbol{\tau}$.
Here $p$ is the hydrostatic pressure.
-We make use of the pressure response throughout this tut (although we refer to it as the pressure).
+We make use of the pressure response throughout this tut (although we refer to it as the pressure).
<h4> Neo-Hookean materials </h4>
One should note that the definitions of the volumetric Kirchhoff stress in the three field formulation
$\boldsymbol{\tau}_{\textrm{vol}} \equiv \widetilde{p} J \mathbf{I}$
- and the subsequent volumetric tangent differs slightly from the general form given in the section on hyperelastic materials where
+ and the subsequent volumetric tangent differs slightly from the general form given in the section on hyperelastic materials where
$\boldsymbol{\tau}_{\textrm{vol}} \equiv p J\mathbf{I}$.
This is because the pressure $\widetilde{p}$ is now a primary field as opposed to a constitutively derived quantity.
-One needs to carefully distinguish between the primary fields and those obtained from the constitutive relations.
+One needs to carefully distinguish between the primary fields and those obtained from the constitutive relations.
@note Although the variables are all expressed in terms of spatial quantities, the domain of integration is the initial configuration.
This approach is called a <em> total-Lagrangian formulation </em>.
-The approach given in step-18, where the domain of integration is the current configuration, could be called an <em> updated Lagrangian formulation </em>.
+The approach given in step-18, where the domain of integration is the current configuration, could be called an <em> updated Lagrangian formulation </em>.
The various merits of these two approaches are discussed widely in the literature.
It should be noted however that they are equivalent.
\int_{\Omega} [-\mathrm{grad}\delta \mathbf{u}:\boldsymbol{\sigma} + \delta \mathbf{u} \cdot\mathbf{b}^\text{p}]~\mathrm{d}v
+ \int_{\partial \Omega} \delta \mathbf{u} \cdot \mathbf{t}^\text{p}~\mathrm{d}a \\
&=
-- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V
+- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V
+ \int_{\Omega_0} \delta \mathbf{u} \cdot J\mathbf{b}^\text{p}~\mathrm{d}V
+ \int_{\partial \Omega_0} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \\
&=
-- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V
+- \int_{\Omega_0} \mathrm{grad}\delta \mathbf{u}:\boldsymbol{\tau}~\mathrm{d}V
+ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V
+ \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \\
&=
-- \int_{\Omega_0} [\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} :\boldsymbol{\tau}~\mathrm{d}V
+- \int_{\Omega_0} [\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} :\boldsymbol{\tau}~\mathrm{d}V
+ \int_{\Omega_0} \delta \mathbf{u} \cdot \mathbf{B}^\text{p}~\mathrm{d}V
+ \int_{\partial \Omega_{0,\sigma}} \delta \mathbf{u} \cdot \mathbf{T}^\text{p}~\mathrm{d}A \, ,
@f}
-where
+where
$[\mathrm{grad}\delta\mathbf{u}]^{\text{sym}} = 1/2[ \mathrm{grad}\delta\mathbf{u} + [\mathrm{grad}\delta\mathbf{u}]^T] $.
We will use an iterative Newton-Raphson method to solve the nonlinear residual equation $R$.
For the sake of simplicity we assume dead loading, i.e. the loading does not change due to the deformation.
The change in a quantity between the known state at $t_{\textrm{n}-1}$
-and the currently unknown state at $t_{\textrm{n}}$ is denoted
+and the currently unknown state at $t_{\textrm{n}}$ is denoted
$\varDelta \{ \bullet \} = { \{ \bullet \} }^{\textrm{n}} - { \{ \bullet \} }^{\textrm{n-1}}$.
-The value of a quantity at the current iteration $\textrm{i}$ is denoted
+The value of a quantity at the current iteration $\textrm{i}$ is denoted
${ \{ \bullet \} }^{\textrm{n}}_{\textrm{i}} = { \{ \bullet \} }_{\textrm{i}}$.
The incremental change between iterations $\textrm{i}$ and $\textrm{i}+1$ is denoted
$d \{ \bullet \} \dealcoloneq \{ \bullet \}_{\textrm{i}+1} - \{ \bullet \}_{\textrm{i}}$.
that is where $\nu \rightarrow 0.5$ (where $\nu$ is <a
href="http://en.wikipedia.org/wiki/Poisson's_ratio">Poisson's ratio</a>), subject to a good choice of the interpolation fields
for $\mathbf{u},~\widetilde{p}$ and $\widetilde{J}$.
-Typically a choice of $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$ is made.
+Typically a choice of $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$ is made.
Here $DGPM$ is the FE_DGPMonomial class.
A popular choice is $Q_1 \times DGPM_0 \times DGPM_0$ which is known as the mean dilatation method (see Hughes (2000) for an intuitive discussion).
This code can accommodate a $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$ formulation.
\end{bmatrix}}_{ \mathbf{\mathsf{F}}(\mathbf{\Xi}_{\textrm{i}}) } \, .
@f}
-There are no derivatives of the pressure and dilatation (primary) variables present in the formulation.
+There are no derivatives of the pressure and dilatation (primary) variables present in the formulation.
Thus the discontinuous finite element interpolation of the pressure and dilatation yields a block
diagonal matrix for
$\mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}$,
-$\mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}$ and
+$\mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{p}}$ and
$\mathbf{\mathsf{K}}_{\widetilde{J}\widetilde{J}}$.
Therefore we can easily express the fields $\widetilde{p}$ and $\widetilde{J}$ on each cell simply
by inverting a local matrix and multiplying it by the local right hand
Note that due to the choice of $\widetilde{p}$ and $\widetilde{J}$ as discontinuous at the element level, all matrices that need to be inverted are defined at the element level.
The procedure to construct the various contributions is as follows:
-- Construct $\mathbf{\mathsf{K}}$.
+- Construct $\mathbf{\mathsf{K}}$.
- Form $\mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}^{-1}$ for element and store where $\mathbf{\mathsf{K}}_{\widetilde{p}\widetilde{J}}$ was stored in $\mathbf{\mathsf{K}}$.
- Form $\overline{\overline{\mathbf{\mathsf{K}}}}$ and add to $\mathbf{\mathsf{K}}_{uu}$ to get $\mathbf{\mathsf{K}}_{\textrm{con}}$
-- The modified system matrix is called ${\mathbf{\mathsf{K}}}_{\textrm{store}}$.
+- The modified system matrix is called ${\mathbf{\mathsf{K}}}_{\textrm{store}}$.
That is
@f[
\mathbf{\mathsf{K}}_{\textrm{store}}
<h3> The material class </h3>
-A good object-oriented design of a Material class would facilitate the extension of this tutorial to a wide range of material types.
+A good object-oriented design of a Material class would facilitate the extension of this tutorial to a wide range of material types.
In this tutorial we simply have one Material class named Material_Compressible_Neo_Hook_Three_Field.
Ideally this class would derive from a class HyperelasticMaterial which would derive from the base class Material.
-The three-field nature of the formulation used here also complicates the matter.
+The three-field nature of the formulation used here also complicates the matter.
The Helmholtz free energy function for the three field formulation is $\Psi = \Psi_\text{vol}(\widetilde{J}) + \Psi_\text{iso}(\overline{\mathbf{b}})$.
The isochoric part of the Kirchhoff stress ${\boldsymbol{\tau}}_{\text{iso}}(\overline{\mathbf{b}})$ is identical to that obtained using a one-field formulation for a hyperelastic material.
Thus, for a three field formulation the constitutive response for the volumetric part of the Kirchhoff stress ${\boldsymbol{\tau}}_{\text{vol}}$ (and the tangent) is not given by the hyperelastic constitutive law as in a one-field formulation.
One can label the term
$\boldsymbol{\tau}_{\textrm{vol}} \equiv \widetilde{p} J \mathbf{I}$
-as the volumetric Kirchhoff stress, but the pressure $\widetilde{p}$ is not derived from the free energy; it is a primary field.
+as the volumetric Kirchhoff stress, but the pressure $\widetilde{p}$ is not derived from the free energy; it is a primary field.
-In order to have a flexible approach, it was decided that the Material_Compressible_Neo_Hook_Three_Field would still be able to calculate and return a volumetric Kirchhoff stress and tangent.
-In order to do this, we choose to store the interpolated primary fields $\widetilde{p}$ and $\widetilde{J}$ in the Material_Compressible_Neo_Hook_Three_Field class associated with the quadrature point.
+In order to have a flexible approach, it was decided that the Material_Compressible_Neo_Hook_Three_Field would still be able to calculate and return a volumetric Kirchhoff stress and tangent.
+In order to do this, we choose to store the interpolated primary fields $\widetilde{p}$ and $\widetilde{J}$ in the Material_Compressible_Neo_Hook_Three_Field class associated with the quadrature point.
This decision should be revisited at a later stage when the tutorial is extended to account for other materials.
That is, the response would be overly stiff.
The initial and final configurations are shown in the image above.
Using symmetry, we solve for only one quarter of the geometry (i.e. a cube with dimension $0.001$).
-The inner-quarter of the upper surface of the domain is subject to a load of $p_0$.
-
-
+The inner-quarter of the upper surface of the domain is subject to a load of $p_0$.
+
+
(detection and stress calculations) itself. An alternative to
additional penalty terms in the free-energy functional would be to
use active set methods such as the one used in step-41.
-- The complete condensation procedure using LinearOperators has been
- coded into the linear solver routine. This could also have been
- achieved through the application of the schur_complement()
+- The complete condensation procedure using LinearOperators has been
+ coded into the linear solver routine. This could also have been
+ achieved through the application of the schur_complement()
operator to condense out one or more of the fields in a more
automated manner.
- Finally, adaptive mesh refinement, as demonstrated in step-6 and
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
Due to the variable coefficient (the curvature there is reduced by the
same factor by which the coefficient is increased), the top region of
the solution is flattened. The gradient of the solution is
-discontinuous along the interface, although this is not very clearly
+discontinuous along the interface, although this is not very clearly
visible in the pictures above. We will look at this in more detail in
the next example.
<h4> Reducing the size of the linear system </h4>
To alleviate the computational cost of solving such large linear systems,
the hybridizable discontinuous Galerkin (HDG) methodology was introduced
-by Cockburn and co-workers (see the references in the recent HDG overview
+by Cockburn and co-workers (see the references in the recent HDG overview
article by Nguyen and Peraire @cite Ngu2012).
The HDG method achieves this goal by formulating the mathematical problem using
skeleton of the triangulation. The global problem is then to solve for the
trace values, which are the only globally coupled unknowns.
<li> Use the known skeleton values as Dirichlet data for solving local
-element-level solutions. This is known as the 'local solver', and is an
+element-level solutions. This is known as the 'local solver', and is an
<i>embarrassingly parallel</i> element-by-element solution process.
</ol>
The above procedure also has a linear algebra interpretation---referred to
as <i>static condensation</i>---that was exploited to reduce the size of the
global linear system by Guyan in the context of continuous Finite Elements
-@cite G65, and by Fraeijs de Veubeke for mixed methods @cite F65. In the
+@cite G65, and by Fraeijs de Veubeke for mixed methods @cite F65. In the
latter case (mixed formulation), the system reduction was achieved through the
use of discontinuous fluxes combined with the introduction of an additional
-auxiliary <i>hybrid</i> variable that approximates the trace of the unknown
+auxiliary <i>hybrid</i> variable that approximates the trace of the unknown
at the boundary of every element. This procedure became known as hybridization
and---by analogy---is the reason why the local discontinuous Galerkin method
-introduced by Cockburn, Gopalakrishnan, and Lazarov in 2009 @cite CGL2009, and
-subsequently developed by their collaborators, eventually came to be known as
+introduced by Cockburn, Gopalakrishnan, and Lazarov in 2009 @cite CGL2009, and
+subsequently developed by their collaborators, eventually came to be known as
the <i>hybridizable discontinuous Galerkin</i> (HDG) method.
-Let us write the complete linear system associated to the HDG problem as a
+Let us write the complete linear system associated to the HDG problem as a
block system with the discrete DG (cell interior) variables $U$ as first block
and the skeleton (face) variables $\Lambda$ as the second block:
@f{eqnarray*}
problem. It solves a small variation of the heat equation discussed first in
step-26 but, since the purpose of this program is only to demonstrate using
more advanced ways to interface with deal.II's time stepping algorithms, only
-solves a simple problem on a uniformly refined mesh.
+solves a simple problem on a uniformly refined mesh.
<h3>Problem statement</h3>
@f}
augmented by appropriate boundary conditions. Here, $v$ is the velocity of
neutrons (for simplicity we assume it is equal to 1 which can be achieved by
-simply scaling the time variable), $D$ is the diffusion coefficient,
+simply scaling the time variable), $D$ is the diffusion coefficient,
$\Sigma_a$ is the absorption cross section, and $S$ is a source. Because we are
only interested in the time dependence, we assume that $D$ and $\Sigma_a$ are
-constant.
+constant.
Since this program only intends to demonstrate how to use advanced time
stepping algorithms, we will only look for the solutions of relatively simple
u_h(x,t) = \sum_j U_j(t) \varphi_j(x)
@f}
with spatial shape functions $\varphi_j(x)$, then multiplying an equation of
-the form
+the form
@f{eqnarray*}
\frac{\partial u(x,t)}{\partial t} = q(t,u(x,t))
@f}
of $q(t,u(x,t))$ (where $q$ is typically the place where spatial
derivatives appear, but this is not of much concern for the moment given that
we only consider time derivatives). In other words, this form fits the general
-scheme above if we write
+scheme above if we write
@f{eqnarray*}
\frac{dy}{dt} = g(t,y) = M^{-1}f(t,y).
@f}
<li> Implicit Runge-Kutta; this class of methods require the solution of a
possibly nonlinear system the stages $k_i$ above, i.e., they have
$a_{ii}\neq 0$ for at least one of the stages $i=1,\ldots,s$.
-</ol>
+</ol>
Many well known time stepping schemes that one does not typically associate
with the names Runge or Kutta can in fact be written in a way so that they,
too, can be expressed in these categories. They oftentimes represent the
lowest-order members of these families.
-<h4>Explicit Runge-Kutta methods</h4>
+<h4>Explicit Runge-Kutta methods</h4>
These methods, only require a function to evaluate $M^{-1}f(t,y)$ but not
(as implicit methods) to solve an equation that involves
Matlab and often abbreviated as RK45 to indicate that the lower and higher order methods
used here are 4th and 5th order Runge-Kutta methods, respectively), Fehlberg,
and Cash-Karp.
-
+
At the time of the writing, only embedded explicit methods have been implemented.
for $y$ in each (sub-)timestep. Internally, this is
done using a Newton-type method and, consequently, they require that the user
provide functions that can evaluate $M^{-1}f(t,y)$ and
-$\left(I-\tau M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently
-$\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M$.
+$\left(I-\tau M^{-1} \frac{\partial f}{\partial y}\right)^{-1}$ or equivalently
+$\left(M - \tau \frac{\partial f}{\partial y}\right)^{-1} M$.
The particular form of this operator results from the fact that each Newton
step requires the solution of an equation of the form
\left(M - \tau \frac{\partial f}{\partial y}\right) \Delta y
= -M h(t,y)
@f}
-for some (given) $h(t,y)$. Implicit methods are
+for some (given) $h(t,y)$. Implicit methods are
always stable, regardless of the time step size, but too large time steps of
course affect the <i>accuracy</i> of the solution, even if the numerical
solution remains stable and bounded.
<ol>
<li> You will need to project the solution to the new mesh when the mesh is changed. Of course,
- the mesh
+ the mesh
used should be the same from the beginning to the end of each time step,
a question that arises because Runge-Kutta methods use multiple
evaluations of the equations within each time step.
One could cite many more examples of complicated geometries in which one wants
to pose and solve a partial differential equation. What this shows is that the
"real" world is much more complicated than what we have shown in almost all of
-the tutorial programs preceding this one.
+the tutorial programs preceding this one.
This program is therefore devoted to showing how one deals with complex
geometries using a concrete application. In particular, what it shows is how
can often be described as a collection of charts (e.g., the northern
hemisphere and the southern hemisphere are each charts) and the domain can then
be describe by an <a
-href="http://en.wikipedia.org/wiki/Atlas_%28topology%29">atlas</a>.
+href="http://en.wikipedia.org/wiki/Atlas_%28topology%29">atlas</a>.
If a domain can be decomposed into an atlas, all we need to do is provide the
pull-back and push-forward functions for each of the charts. In deal.II, this
<img src="https://www.dealii.org/images/steps/developer/step-53.topo.png" alt="">
In particular, though you cannot see this here, the top surface is not
-just colored by the elevation but is, in fact, deformed to follow the
+just colored by the elevation but is, in fact, deformed to follow the
correct topography.
While the actual application is not relevant here, the geometry is. The domain
we are interested in is a part of the Earth that ranges from the surface to a
depth of 500km, from 26 to 35 degrees East of the Greenwich meridian, and from
5 degrees North of the equator to 10 degrees South.
-This description of the geometry suggests to start with a box
+This description of the geometry suggests to start with a box
$\hat U=[26,35]\times[-10,5]\times[-500000,0]$ (measured in degrees,
degrees, and meters) and to provide a map $\varphi$ so
that $\varphi^{-1}(\hat U)=\Omega$ where $\Omega$ is the domain we
this using the <a href="http://en.wikipedia.org/wiki/WGS84">WGS 84</a>
reference standard for the Earth shape. The formula used in WGS 84 to obtain
a position in Cartesian coordinates from longitude, latitude, and elevation
- is
+ is
@f[
\mathbf x = \varphi_\text{WGS84}^{-1}(\phi,\theta,d)
= \left(
\right),
@f]
where $\bar R(\theta)=\frac{R}{\sqrt{1-(e \sin\theta)^2}}$, and radius and
- ellipticity are given by $R=6378137\text{m}, e=0.081819190842622$. In this formula,
+ ellipticity are given by $R=6378137\text{m}, e=0.081819190842622$. In this formula,
we assume that the arguments to sines and cosines are evaluated in degree, not
radians (though we will have to change this assumption in the code).
this using real topography data (see below for a description of where
this data comes from). Using this data set, we can look up elevations on a
latitude-longitude mesh laid over the surface of the Earth. Starting with
- the box $\hat U=[26,35]\times[-10,5]\times[-500000,0]$, we will therefore
+ the box $\hat U=[26,35]\times[-10,5]\times[-500000,0]$, we will therefore
first stretch it in vertical direction before handing it off to the WGS 84
function: if $h(\hat\phi,\hat\theta)$ is the height at longitude $\hat\phi$
and latitude $\hat\theta$, then we define
Using these two functions, we can then define the entire push-forward function
$\varphi^{-1}: \hat U \rightarrow \Omega$ as
@f[
- \mathbf x
- =
+ \mathbf x
+ =
\varphi^{-1}(\hat\phi,\hat\theta,\hat d)
=
\varphi_\text{WGS84}^{-1}(\varphi_\text{topo}^{-1}(\hat\phi,\hat\theta,\hat d)).
pull-back operation, which we can write as
@f[
(\hat\phi,\hat\theta,\hat d)
- =
+ =
\varphi(\mathbf x)
=
\varphi_\text{topo}(\varphi_\text{WGS84}(\mathbf x)).
<h2>Implementation</h2>
-There are a number of issues we need to address in the program. At the largest scale,
+There are a number of issues we need to address in the program. At the largest scale,
we need to write a class that implements the interface of ChartManifold. This involves
a function <code>push_forward()</code> that takes a point
in the reference domain $\hat U$ and transform it into real space using the function
@endcode
The transformations above have two parts: the WGS 84 transformations and the topography
-transformation. Consequently, the <code>AfricaGeometry</code> class will have
-additional (non-virtual) member functions
+transformation. Consequently, the <code>AfricaGeometry</code> class will have
+additional (non-virtual) member functions
<code>AfricaGeometry::push_forward_wgs84()</code> and
<code>AfricaGeometry::push_forward_topo()</code> that implement these two pieces, and
corresponding pull back functions.
@endcode
The data is formatted as <code>latitude longitude elevation</code> where the first two
columns are provided in degrees North of the equator and degrees East of the Greenwich
-meridian. The final column is given in meters above the WGS 84 zero elevation.
+meridian. The final column is given in meters above the WGS 84 zero elevation.
In the transformation functions, we need to evaluate $h(\hat\phi,\hat\theta)$ for a given
longitude $\hat\phi$ and latitude $\hat\theta$. In general, this data point will not be
prone. Fortunately, we can somehow shoehorn this data set into an existing class:
Functions::InterpolatedUniformGridData . Unfortunately, the class does not fit the bill
quite exactly and so we need to work around it a bit. The problem comes from the way
-we initialize this class: in its simplest form, it takes a stream of values that it
+we initialize this class: in its simplest form, it takes a stream of values that it
assumes form an equispaced mesh in the $x-y$ plane (or, here, the $\phi-\theta$ plane).
Which is what they do here, sort of: they are ordered latitude first, longitude second;
and more awkwardly, the first column starts at the largest values and counts down,
-rather than the usual other way around.
+rather than the usual other way around.
Now, while tutorial programs are meant to illustrate how to code with deal.II, they do
not necessarily have to satisfy the same quality standards as one would have to do
On the other hand, tutorial programs are best if they are short and demonstrate key
points rather than dwell on unimportant aspects and, thereby, obscure what we really
want to show. Consequently, we will allow ourselves a bit of leeway:
-- since this program is intended solely for a particular geometry around the area
- of the East-African rift and since this is precisely the area described by the data
- file, we will hardcode in the program that there are
+- since this program is intended solely for a particular geometry around the area
+ of the East-African rift and since this is precisely the area described by the data
+ file, we will hardcode in the program that there are
$1139\times 660$ pieces of data;
-- we will hardcode the boundaries of the data
+- we will hardcode the boundaries of the data
$[-11.98333^\circ,6.983333^\circ]\times[25^\circ,35.98333^\circ]$;
- we will lie to the Functions::InterpolatedUniformGridData class: the class will
only see the data in the last column of this data file, and we will pretend that
interpolating table class for a value at $-\hat\theta$. With this little
trick, we can avoid having to switch around the order of data as read from
file.
-
+
All of this then calls for a class that essentially looks like this:
@code
class AfricaTopography
sphere (though if one modified the program so that it does produce a sphere and
looked at them at the same time, the difference between the overall sphere and
WGS 84 shape is quite apparent). Apparently, Earth is actually quite a flat place.
-Of course we already know this from satellite pictures.
+Of course we already know this from satellite pictures.
However, we can tease out something more by
coloring cells by their volume. This both produces slight variations in hue
along the top surface and something for the visualization programs to apply
A zoom-in of this picture shows the vertical displacement quite clearly (here,
looking from the West-Northwest over the rift valley, the triple peaks
-of
+of
<a href="http://en.wikipedia.org/wiki/Mount_Stanley">Mount Stanley</a>,
<a href="http://en.wikipedia.org/wiki/Mount_Speke">Mount Speke</a>, and
<a href="http://en.wikipedia.org/wiki/Mount_Baker_%28Uganda%29">Mount Baker</a>
<img src="https://www.dealii.org/images/steps/developer/step-53.topozoom.png" alt="">
-These image were produced with three small modifications:
+These image were produced with three small modifications:
<ol>
<li> An additional seventh mesh refinement towards the top surface for the
first of these two pictures, and a total of nine for the second. In the
under 1km in vertical direction. (The picture was also created using a
more resolved data set; however, it is too big to distribute as part of
the tutorial.)
-
- <li> The addition of the following function that, given a point
+
+ <li> The addition of the following function that, given a point
<code>x</code> computes the elevation by converting the point to
reference WGS 84 coordinates and only keeping the depth variable (the
function is, consequently, a simplified version of the
5,
ScalarFunctionFromFunctionObject<3>(get_elevation),
boundary_values);
- for (std::map<unsigned int,double>::const_iterator p = boundary_values.begin();
+ for (std::map<unsigned int,double>::const_iterator p = boundary_values.begin();
p!=boundary_values.end(); ++p)
elevation[p->first] = p->second;
}
like an object of a class derived from the Function class that we want
to use in VectorTools::interpolate_boundary_values(). Having so gotten a list
of degrees of freedom located at the top boundary and corresponding elevation
-values, we just go down this list and set these elevations in the
+values, we just go down this list and set these elevations in the
<code>elevation</code> vector (leaving all interior degrees of freedom at
their original zero value). This vector is then output using DataOut as
usual and can be visualized as shown above.
to access all the sub-shapes (such as vertices, edges and faces) composing it,
along with their geometrical description. In the deal.II framework, the
topological entities composing a shape are used to create a corresponding
-Manifold representation. In step-6 we saw how to use GridGenerator::hyper_sphere()
+Manifold representation. In step-6 we saw how to use GridGenerator::hyper_sphere()
to create a hyper sphere, which automatically attaches a SphericalManifold
to all boundary faces. This guarantees that boundary faces stay on a
sphere or circle during mesh refinement. The functions of the CAD modeling interface
<h1>Results</h1>
-The program execution produces a series of mesh files <code>3d_mesh_*.vtk</code>
+The program execution produces a series of mesh files <code>3d_mesh_*.vtk</code>
that we can visualize with any of the usual visualization programs that can read the VTK
-file format.
+file format.
The following table illustrates the results obtained employing the normal projection strategy. The first two
rows of the table show side views of the grids obtained for progressive levels
<table style="width:90%">
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.common_0.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_1.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_1.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_2.png" alt="" width="400"></td>
</tr>
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_3.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_4.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_4.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_5.png" alt="" width="400"></td>
</tr>
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_front_3.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_front_4.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_front_4.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_front_5.png" alt="" width="400"></td>
</tr>
</table>
<table style="width:90%">
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.common_0.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_1.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_1.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_2.png" alt="" width="400"></td>
</tr>
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_3.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_4.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_4.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_5.png" alt="" width="400"></td>
</tr>
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_front_3.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_front_4.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_front_4.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.directional_front_5.png" alt="" width="400"></td>
</tr>
</table>
The images confirm that the quality of the mesh obtained with a directional projection is sensibly higher than that obtained projecting along the
surface normal. Yet, a number of elements elongated in the y-direction are observed around the bottom of the bulb, where the surface is almost parallel to the
-direction chosen for the projection.
+direction chosen for the projection.
The final test shows results using instead the projection normal to the faces:
<table style="width:90%">
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.common_0.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_1.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_1.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_2.png" alt="" width="400"></td>
</tr>
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_3.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_4.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_4.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_5.png" alt="" width="400"></td>
</tr>
<tr>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_3.png" alt="" width="400"></td>
- <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_4.png" alt="" width="400"></td>
+ <td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_4.png" alt="" width="400"></td>
<td><img src="https://www.dealii.org/images/steps/developer/step-54.normal_to_mesh_front_5.png" alt="" width="400"></td>
</tr>
</table>
The pictures confirm that the normal to mesh projection approach leads to grids that remain evenly spaced
throughtout the refinement steps. At the same time, these meshes represent rather well the original geometry even in the bottom region
of the bulb, which is not well recovered employing the directional projector or the normal projector.
-
+
@dealiiTutorialDOI{10.5281/zenodo.400995,https://zenodo.org/badge/DOI/10.5281/zenodo.400995.svg}
-<a name="Intro"></a>
+<a name="Intro"></a>
<h1>Introduction</h1>
<h3> Stokes Problem </h3>
@f{eqnarray*}
- 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u})
+ (\nabla \textbf{u})^T\right] + \nabla p & =& f \\
- - \nabla \cdot u &=& 0
+ - \nabla \cdot u &=& 0
@f}
Note that we are using the deformation tensor instead of $\Delta u$ (a
The weak form of
the discrete equations naturally leads to the following linear system
for the nodal values of the velocity and pressure fields:
-@f{eqnarray*}
+@f{eqnarray*}
\left(\begin{array}{cc} A & B^T \\ B & 0
\end{array}\right) \left(\begin{array}{c} U \\ P \end{array}\right) =
-\left(\begin{array}{c} F \\ 0 \end{array}\right).
+\left(\begin{array}{c} F \\ 0 \end{array}\right).
@f}
Our goal is to compare several solution approaches. While step-22
block system at once using FMGRES with an efficient
preconditioner, in the spirit of the approach outlined in the "Results"
section of step-22. The idea is as follows: if we find a block
-preconditioner $P$ such that the matrix
+preconditioner $P$ such that the matrix
@f{eqnarray*}
-\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1}
+\left(\begin{array}{cc} A & B^T \\ B & 0 \end{array}\right) P^{-1}
@f}
is simple, then an iterative solver with that preconditioner will
@f{eqnarray*}
P^{-1} = \left(\begin{array}{cc} A & B^T \\ 0 &
- S \end{array}\right)^{-1}
+ S \end{array}\right)^{-1}
@f}
is a good choice. Let $\widetilde{A^{-1}}$ be an approximation of $A^{-1}$
and $\widetilde{S^{-1}}$ of $S^{-1}$, we see
-@f{eqnarray*}
+@f{eqnarray*}
P^{-1} =
\left(\begin{array}{cc} A^{-1} & 0 \\ 0 & I \end{array}\right)
\left(\begin{array}{cc} I & B^T \\ 0 & -I \end{array}\right)
of GMG.
As a comparison, instead of FGMRES, we also use the direct solver
-UMFPACK on the whole system to compare our results with. If you want to use
+UMFPACK on the whole system to compare our results with. If you want to use
a direct solver (like UMFPACK), the system needs to be invertible. To avoid
the one dimensional null space given by the constant pressures, we fix the first pressure unknown
to zero. This is not necessary for the iterative solvers.
f$ such that:
@f{align*}
-{\bf f} = - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p.
+{\bf f} = - 2 \text{div} \frac {1}{2} \left[ (\nabla \textbf{u}) + (\nabla \textbf{u})^T\right] + \nabla p.
@f}
Using the reference solution above, we obtain:
-@f{eqnarray*}
+@f{eqnarray*}
{\bf f} &=& (2 \pi^2 \sin (\pi x),- \pi^3 y \cos(\pi
x),- \pi^3 z \cos(\pi x))\\ & & + (\pi \cos(\pi x) \cos(\pi y)
\sin(\pi z) ,- \pi \sin(\pi y) \sin(\pi x) \sin(\pi z), \pi \cos(\pi
z) \sin(\pi x) \cos(\pi y)) @f}
-<h3> Computing Errors </h3>
+<h3> Computing Errors </h3>
Because we do not enforce the mean
pressure to be zero for our numerical solution in the linear system,
<h3> Differences from the Step 22 tutorial </h3>
-
+
The main difference between step-56 and step-22 is that we use block
solvers instead of the Schur Complement approach used in
step-22. Details of this approach can be found under the "Block Schur
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
used to be implemented in FORTRAN, a language that does not care about
whether something is spelled in lower or UPPER case letters, with programmers
often choosing upper case letters habitually), the loop above is often
-referenced in publications about mesh adaptivity as the
+referenced in publications about mesh adaptivity as the
SOLVE-ESTIMATE-MARK-REFINE loop (with this spelling).
Beyond this structure, however, there are a variety of ways to achieve
neighbor to the opposite vertex; otherwise a slightly more involved
construction is necessary that adds more new vertices on at least one
other edge of the neighboring cell, and then may propagate to the neighbors
- of the neighbor until the algorithm terminates. This is hard to describe
+ of the neighbor until the algorithm terminates. This is hard to describe
in words, and because deal.II does not use triangles not worth the time here.
But if you're curious, you can always watch video lecture 15 at the link
shown at the top of this introduction.
to elliptic only problems. Although it will create non-optimal meshes for other
equations, it is often a good way to quickly produce meshes that are
well adapted to the features of solutions, such as regions of great
-variation or discontinuities.
+variation or discontinuities.
The AffineConstraints class can handle such constraints as well, which makes it
convenient to let the same object we use for hanging node constraints
also deal with these Dirichlet boundary conditions.
-This way, we don't need to apply the boundary conditions after assembly
+This way, we don't need to apply the boundary conditions after assembly
(like we did in the earlier steps).
All that is necessary is that we call the variant of
VectorTools::interpolate_boundary_values() that returns its information
in an AffineConstraints object, rather than the `std::map` we have used
in previous tutorial programs.
-
+
<h3> Other things this program shows </h3>
-
-
+
+
Since the concepts used for locally refined grids are so important,
we do not show much other material in this example. The most
important exception is that we show how to use biquadratic elements
latter kind would consider the gradient (the "Jacobian") of the mapping
from the reference cell $\hat K=[0,1]^d$ to the real cell $K$; this
gradient is a matrix, and a quantity that enters into error estimates
-is the maximum over all points on the reference cell of the ratio of the
+is the maximum over all points on the reference cell of the ratio of the
largest to the smallest eigenvalue of this matrix. It is again not difficult
to see that this ratio is constant if the cell $K$ is an affine image of
$\hat K$, and that it is one for squares and cubes.
still in the <a href="https://en.wikipedia.org/wiki/Sobolev_space">Sobolev space</a>
$W^{1,\infty}$ (i.e., roughly speaking, in the
space of functions whose derivatives are bounded). On the other hand,
-we know that in the most
+we know that in the most
extreme cases -- i.e., where the domain has reentrant corners, the
right hand side only satisfies $f\in H^{-1}$, or the coefficient $a$ is only in
$L^\infty$ -- all we can expect is that $u\in H^1$ (i.e., the
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
<i>
This program was contributed by Zhuoran Wang.
Some more information about this program, as well as more numerical
-results, are presented in @cite Wang2019 .
+results, are presented in @cite Wang2019 .
</i>
<a name="Intro"></a>
<h1>Introduction</h1>
-This tutorial program presents an implementation of the "weak Galerkin"
+This tutorial program presents an implementation of the "weak Galerkin"
finite element method for the Poisson equation. In some sense, the motivation for
considering this method starts from the same point as in step-51: We would like to
consider discontinuous shape functions, but then need to address the fact that
to rewriting the equation as a first-order system, but it is also not quite
as easy to implement. However, as we will see in the following, this
additional effort is not prohibitive.
-
+
<h3> Weak Galerkin finite element methods </h3>
-Weak Galerkin Finite Element Methods (WGFEMs) use discrete weak functions
-to approximate scalar unknowns, and discrete weak gradients to
-approximate classical gradients.
-The method was originally introduced by Junping Wang and Xiu Ye
-in the paper
+Weak Galerkin Finite Element Methods (WGFEMs) use discrete weak functions
+to approximate scalar unknowns, and discrete weak gradients to
+approximate classical gradients.
+The method was originally introduced by Junping Wang and Xiu Ye
+in the paper
<a href="https://doi.org/10.1016/j.cam.2012.10.003">
-<i>A weak Galerkin finite element method for second order elliptic problems</i>,
-J. Comput. Appl. Math., 103-115, 2013</a>.
-Compared to the continuous Galerkin method,
-the weak Galerkin method satisfies important physical properties, namely
-local mass conservation and bulk normal flux continuity.
-It results in a SPD linear system, and optimal convergence rates can
+<i>A weak Galerkin finite element method for second order elliptic problems</i>,
+J. Comput. Appl. Math., 103-115, 2013</a>.
+Compared to the continuous Galerkin method,
+the weak Galerkin method satisfies important physical properties, namely
+local mass conservation and bulk normal flux continuity.
+It results in a SPD linear system, and optimal convergence rates can
be obtained with mesh refinement.
<h3> The equation to solve </h3>
-This program solves the Poisson equation
+This program solves the Poisson equation
using the weak Galerkin finite element method:
@f{align*}{
\nabla \cdot \left( -\mathbf{K} \nabla p \right)
\qquad \mathbf{x} \in \Omega, \\
p &= p_D,\qquad \mathbf{x} \in \Gamma^D, \\
\mathbf{u} \cdot \mathbf{n} &= u_N,
- \qquad \mathbf{x} \in \Gamma^N,
+ \qquad \mathbf{x} \in \Gamma^N,
@f}
-where $\Omega \subset \mathbb{R}^n (n=2,3)$ is a bounded domain.
+where $\Omega \subset \mathbb{R}^n (n=2,3)$ is a bounded domain.
In the context of the flow of a fluid through a porous medium,
-$p$ is the pressure, $\mathbf{K}$ is a permeability tensor,
+$p$ is the pressure, $\mathbf{K}$ is a permeability tensor,
$f$ is the source term, and
$p_D, u_N$ represent Dirichlet and Neumann boundary conditions.
-We can introduce a flux, $\mathbf{u} = -\mathbf{K} \nabla p$, that corresponds
+We can introduce a flux, $\mathbf{u} = -\mathbf{K} \nabla p$, that corresponds
to the Darcy velocity (in the way we did in step-20) and this variable will
be important in the considerations below.
-In this program, we will consider a test case where the exact pressure
+In this program, we will consider a test case where the exact pressure
is $p = \sin \left( \pi x\right)\sin\left(\pi y \right)$ on the unit square domain,
with homogeneous Dirichelet boundary conditions and $\mathbf{K}$ the identity matrix.
Then we will calculate $L_2$ errors of pressure, velocity, and flux.
\mathcal{A}\left(p,q\right)
\dealcoloneq \int_\Omega \left(\mathbf{K} \nabla p\right) \cdot \nabla q \;\mathrm{d}x,
@f}
-and
+and
@f{equation*}
\mathcal{F}\left(q\right)
\dealcoloneq \int_\Omega f \, q \;\mathrm{d}x
- \int_{\Gamma^N} u_N q \; \mathrm{d}x.
@f}
-Here, we have integrated by parts in the bilinear form, and we are evaluating
+Here, we have integrated by parts in the bilinear form, and we are evaluating
the gradient of $p,p$ in the interior and the values of $q$ on the boundary
of the domain. All of this is well defined because we assume that the solution
is in $H^1$ for which taking the gradient and evaluating boundary values
-are valid operations.
+are valid operations.
The idea of the weak Galerkin method is now to approximate the exact $p$
solution with a <i>discontinuous function</i> $p_h$. This function may only be
\dealcoloneq \sum_{K \in \mathbb{T}}
\int_K \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q_h \;\mathrm{d}x,
@f}
-and
+and
@f{equation*}
\mathcal{F}\left(q_h\right)
\dealcoloneq \sum_{K \in \mathbb{T}} \int_K f \, q_h^\circ \;\mathrm{d}x
The key point is that here, we have replaced the gradient $\nabla p_h$ by the
<i>discrete weak gradient</i> operator
$\nabla_{w,d} p_h$ that makes sense for our peculiarly defined approximation $p_h$.
-
+
The question is then how that operator works. For this, let us first say how we
think of the discrete approximation $p_h$ of the pressure. As mentioned above,
the "function" $p_h$ actually consists of two parts: the values $p_h^\circ$ in
think of $\nabla_{w,d} p_h$ of some function defined on $K$ that approximates
the gradient; in particular, if $p_h|_K$ was the restriction of a differentiable
function (to the interior and boundary of $K$ -- which would make it continuous
-between the interior and boundary), then
+between the interior and boundary), then
$\nabla_{w,d} p_h$ would simply be the exact gradient $\nabla p_h$. But, since
$p_h|_K$ is not continuous between interior and boundary of $K$, we need a more
general definition; furthermore, we can not deal with arbitrary functions, and
DGQ_k(K) \times DGQ_r(\partial K) \rightarrow RT_s(K)$ (where $RT_s(K)$ is the
vector-valued Raviart-Thomas space of order $s$ on cell $K$) in the following way:
@f{equation*}{
- \int_K \mathbf v_h \cdot (\nabla_{w,d} p_h)
+ \int_K \mathbf v_h \cdot (\nabla_{w,d} p_h)
=
-\int_K (\nabla \cdot \mathbf v_h) p_h^\circ
+\int_{\partial K} (\mathbf v_h \cdot \mathbf n) p_h^\partial,
-@f}
+@f}
for all test functions $\mathbf v_h \in RT_s(K)$.
This is, in essence, simply an application of the integration-by-parts
formula. In other words, for a given $p_h=(p^\circ_h,p^\partial_h)$,
we need to think of $\nabla_{w,d} p_h|_K$ as that
Raviart-Thomas function of degree $s$ for which the left hand side and right hand side
-are equal for all test functions.
+are equal for all test functions.
@note It may be worth pointing out that while the weak discrete
gradient is an element of the Raviart-Thomas space $RT_s(K)$ on each
as we always do, i.e., we can write
@f{equation*}{
p_h(\mathbf x) = \sum_j P_j \varphi_j(\mathbf x).
-@f}
+@f}
Here, since $p_h$ has two components (the interior and the interface components),
the same must hold true for the basis functions $\varphi_j(\mathbf x)$, which we
can write as $\varphi_j = (\varphi_j^\circ,\varphi_j^\partial)$. If you've
shape functions will be of either one or the other kind. That is not important,
here, however. What is important is that we need to wonder how we can represent
$\nabla_{w,d} \varphi_j$ because that is clearly what will appear in the
-problem when we want to implement the bilinear form
+problem when we want to implement the bilinear form
@f{equation*}
\mathcal{A}_h\left(p_h,q_h\right)
= \sum_{K \in \mathbb{T}}
= \sum_k C_{jk}^K \mathbf v_k|_K
@f}
where the functions $\mathbf v_k \in DGRT_s$, and where $C^K$ is a matrix of
-dimension
+dimension
@f{align*}{
\text{dim}\left(DGQ_k(K) \times DGQ_r(K)\right) &\times \text{dim}\left(RT_s(K)\right)
\\
Using this expansion, we can easily use the definition of the weak
discrete gradient above to define what the matrix is going to be:
@f{equation*}{
- \int_K \mathbf v_i \cdot \left(\sum_k C_{jk}^K \mathbf v_k\right)
+ \int_K \mathbf v_i \cdot \left(\sum_k C_{jk}^K \mathbf v_k\right)
=
-\int_K (\nabla \cdot \mathbf v_i) \varphi_j^\circ
+\int_{\partial K} (\mathbf v_i \cdot \mathbf n) \varphi_j^\partial,
This clearly leads to a linear system of the form
@f{equation*}{
- \sum_k M_{ik}^K C_{jk}^K
+ \sum_k M_{ik}^K C_{jk}^K
=
G_{ij}^K
@f}
with
@f{equation*}{
M_{ik}^K = \int_K \mathbf v_i \cdot \mathbf v_k,
- \qquad\qquad
+ \qquad\qquad
G_{ij}^K = -\int_K (\nabla \cdot \mathbf v_i) \varphi_j^\circ
+\int_{\partial K} (\mathbf v_i \cdot \mathbf n) \varphi_j^\partial,
@f}
-and consequently
+and consequently
@f{equation*}{
\left(C^K\right)^T = \left(M^K\right)^{-1} G^K.
@f}
(In this last step, we have assumed that the indices $i,j,k$ only range
over those degrees of freedom active on cell $K$,
-thereby ensuring that the mass matrix on the space $RT_s(K)$ is invertible.)
+thereby ensuring that the mass matrix on the space $RT_s(K)$ is invertible.)
Equivalently, using the symmetry of the matrix $M$, we have that
@f{equation*}{
C^K = \left(G^K\right)^{T} \left(M^K\right)^{-1}.
@f{equation*}{
A^K_{ij} = \int_K \left({\mathbf K} \nabla_{w,d} \varphi_i\right) \cdot \nabla_{w,d} \varphi_j.
@f}
-As explained above, we can expand $\nabla_{w,d} \varphi_i$ in terms of the
+As explained above, we can expand $\nabla_{w,d} \varphi_i$ in terms of the
Raviart-Thomas basis on each cell, and similarly for $\nabla_{w,d} \varphi_j$:
@f{equation*}{
- A^K_{ij} = \int_K
+ A^K_{ij} = \int_K
\left(
- {\mathbf K}
+ {\mathbf K}
\sum_k C_{ik}^K \mathbf v_k|_K
- \right)
- \cdot
+ \right)
+ \cdot
\sum_l C_{jl}^K \mathbf v_l|_K.
@f}
By re-arranging sums, this yields the following expression:
@f{equation*}{
A^K_{ij} =
\sum_k \sum_l C_{ik}^K C_{jl}^K
- \int_K
+ \int_K
\left(
- {\mathbf K}
+ {\mathbf K}
\mathbf v_k|_K
- \right)
- \cdot
+ \right)
+ \cdot
\mathbf v_l|_K.
@f}
So, if we have the matrix $C^K$ for each cell $K$, then we can easily compute
Here,
@f{equation*}{
H^K_{kl} =
- \int_K
+ \int_K
\left(
- {\mathbf K}
+ {\mathbf K}
\mathbf v_k|_K
- \right)
- \cdot
+ \right)
+ \cdot
\mathbf v_l|_K,
@f}
which is really just the mass matrix on cell $K$ using the Raviart-Thomas
on each element the gradient of the numerical pressure $\nabla p$ can be
approximated by discrete weak gradients $ \nabla_{w,d}\phi_i$:
@f{equation*}
-\nabla_{w,d} p_h
-= \nabla_{w,d} \left(\sum_{i} P_i \phi_i\right)
+\nabla_{w,d} p_h
+= \nabla_{w,d} \left(\sum_{i} P_i \phi_i\right)
= \sum_{i} P_i \nabla_{w,d}\phi_i.
@f}
On cell $K$,
the numerical velocity $ \mathbf{u}_h = -\mathbf{K} \nabla_{w,d}p_h$ can be written as
@f{align*}{
- \mathbf{u}_h
- &= -\mathbf{K} \nabla_{w,d} p_h
+ \mathbf{u}_h
+ &= -\mathbf{K} \nabla_{w,d} p_h
= -\mathbf{K}\sum_{i} \sum_{j} P_i C^K_{ij}\mathbf{v}_j,
@f}
where $C^K$ is the expansion matrix from above, and
$\mathbf{v}_j$ is the basis function of the $RT$ space on a cell.
Unfortunately, $\mathbf{K} \mathbf{v}_j$ may not be in the $RT$ space
-(unless, of course, if $\mathbf K$ is constant times the identity matrix).
+(unless, of course, if $\mathbf K$ is constant times the identity matrix).
So, in order to represent it in a finite element program, we need to
project it back into a finite dimensional space we can work with. Here,
we will use the $L_2$-projection to project it back to the (broken) $RT$
-space.
-
-We define the projection as
-$ \mathbf{Q}_h \left( \mathbf{K}\mathbf{v}_j \right) =
-\sum_{k} d_{jk}\mathbf{v}_k$ on each cell $K$.
-For any $j$,
-$\left( \mathbf{Q}_h \left( \mathbf{Kv}_j \right),\mathbf{v}_k \right)_K =
-\left( \mathbf{Kv}_j,\mathbf{v}_k \right)_K.$
+space.
+
+We define the projection as
+$ \mathbf{Q}_h \left( \mathbf{K}\mathbf{v}_j \right) =
+\sum_{k} d_{jk}\mathbf{v}_k$ on each cell $K$.
+For any $j$,
+$\left( \mathbf{Q}_h \left( \mathbf{Kv}_j \right),\mathbf{v}_k \right)_K =
+\left( \mathbf{Kv}_j,\mathbf{v}_k \right)_K.$
So, rather than the formula shown above, the numerical velocity on cell $K$
instead becomes
@f{equation*}
-\mathbf{u}_h = \mathbf{Q}_h \left( -\mathbf{K}\nabla_{w,d}p_h \right) =
+\mathbf{u}_h = \mathbf{Q}_h \left( -\mathbf{K}\nabla_{w,d}p_h \right) =
-\sum_i \sum_j P_i B^K_{ij}\mathbf{Q}_h \left( \mathbf{K}\mathbf{v}_j \right),
@f}
and we have the following system to solve for the coefficients $d_{jk}$:
=
\left( \mathbf{Kv}_j,\mathbf{v}_k \right).
@f}
-In the implementation below, the matrix with elements
+In the implementation below, the matrix with elements
$
d_{jk}
$
is called <code>cell_matrix_D</code>,
-whereas the matrix with elements
+whereas the matrix with elements
$
\left( \mathbf{Kv}_j,\mathbf{v}_k \right)
$
Then the elementwise velocity is
@f{equation*}
-\mathbf{u}_h = -\sum_{i} \sum_{j}P_ic_{ij}\sum_{k}d_{jk}\mathbf{v}_k =
+\mathbf{u}_h = -\sum_{i} \sum_{j}P_ic_{ij}\sum_{k}d_{jk}\mathbf{v}_k =
\sum_{k}- \left(\sum_{j} \sum_{i} P_ic_{ij}d_{jk} \right)\mathbf{v}_k,
@f}
-where $-\sum_{j} \sum_{i} P_ic_{ij}d_{jk}$ is called
+where $-\sum_{j} \sum_{i} P_ic_{ij}d_{jk}$ is called
`cell_velocity` in the code.
Using this velocity obtained by "postprocessing" the solution, we can
-define the $L_2$-errors of pressure, velocity, and flux
+define the $L_2$-errors of pressure, velocity, and flux
by the following formulas:
@f{align*}{
\|p-p_h^\circ\|^2
&= \sum_{K \in \mathbb{T}} \sum_{\gamma \subset \partial K}
\frac{|K|}{|\gamma|} \|\mathbf{u} \cdot \mathbf{n} - \mathbf{u}_h \cdot \mathbf{n}\|_{L_2(\gamma)}^2,
@f}
-where $| K |$ is the area of the element,
-$\gamma$ are faces of the element,
+where $| K |$ is the area of the element,
+$\gamma$ are faces of the element,
$\mathbf{n}$ are unit normal vectors of each face. The last of these
norms measures the accuracy of the normal component of the velocity
vectors over the interfaces between the cells of the mesh. The scaling
</tr>
</table>
-From the figures, we can see that with the mesh refinement, the maximum and
-minimum pressure values are approaching the values we expect.
-Since the mesh is a rectangular mesh and numbers of cells in each direction is even, we
-have symmetric solutions. From the 3d figures on the right,
+From the figures, we can see that with the mesh refinement, the maximum and
+minimum pressure values are approaching the values we expect.
+Since the mesh is a rectangular mesh and numbers of cells in each direction is even, we
+have symmetric solutions. From the 3d figures on the right,
we can see that on $\mbox{WG}(Q_0,Q_0;RT_{[0]})$, the pressure is a constant
in the interior of the cell, as expected.
<h4>Convergence table</h4>
We run the code with differently refined meshes (chosen in the `make_grid()` function)
-and get the following convergence rates of pressure,
+and get the following convergence rates of pressure,
velocity, and flux (as defined in the introduction).
<table align="center" class="doxtable">
We can repeat the experiment from above using the next higher polynomial
degree:
-The following figures are interior pressures and face pressures implemented using
-$\mbox{WG}(Q_1,Q_1;RT_{[1]})$. The mesh is refined 4 times. Compared to the
+The following figures are interior pressures and face pressures implemented using
+$\mbox{WG}(Q_1,Q_1;RT_{[1]})$. The mesh is refined 4 times. Compared to the
previous figures using
$\mbox{WG}(Q_0,Q_0;RT_{[0]})$, on each cell, the solution is no longer constant
on each cell, as we now use bilinear polynomials to do the approximation.
<h4>Convergence table</h4>
-The following are the convergence rates of pressure, velocity, and flux
+The following are the convergence rates of pressure, velocity, and flux
we obtain from using the $\mbox{WG}(Q_1,Q_1;RT_{[1]})$ element combination:
<table align="center" class="doxtable">
<h3>Test results on <i>WG(Q<sub>2</sub>,Q<sub>2</sub>;RT<sub>[2]</sub>)</i></h3>
Let us go one polynomial degree higher.
-The following are interior pressures and face pressures implemented using
+The following are interior pressures and face pressures implemented using
$WG(Q_2,Q_2;RT_{[2]})$, with mesh size $h = 1/32$ (i.e., 5 global mesh
refinement steps). In the program, we use
`data_out_face.build_patches(fe.degree)` when generating graphical output
<h4>Convergence table</h4>
As before, we can generate convergence data for the
-$L_2$ errors of pressure, velocity, and flux
+$L_2$ errors of pressure, velocity, and flux
using the $\mbox{WG}(Q_2,Q_2;RT_{[2]})$ combination:
<table align="center" class="doxtable">
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
First we create a HDF5 file where we store the parameters and the results of
the simulation.
-Each of the simulations (displacement and calibration) is stored in a separate HDF5 group:
+Each of the simulations (displacement and calibration) is stored in a separate HDF5 group:
@code{.py}
import numpy as np
import h5py
group.attrs['dimension_x'] = 2e-5
# The waveguide width is equal to dimension_y
group.attrs['dimension_y'] = 2e-8
-
+
# Position of the probe that we use to measure the flux
group.attrs['probe_pos_x'] = 8e-6
group.attrs['probe_pos_y'] = 0
group.attrs['probe_width_y'] = 2e-08
-
+
# Number of points in the probe
group.attrs['nb_probe_points'] = 5
-
+
# Global refinement
group.attrs['grid_level'] = 1
if group == displacement:
group.attrs['material_b_rho'] = 2000
else:
- group.attrs['material_b_rho'] = 3200
+ group.attrs['material_b_rho'] = 3200
group.attrs['lambda'] = (group.attrs['youngs_modulus'] * group.attrs['poissons_ratio'] /
((1 + group.attrs['poissons_ratio']) *
(1 - 2 * group.attrs['poissons_ratio'])))
else:
group.attrs['simulation_name'] = 'phononic_cavity_calibration'
group.attrs['save_vtu_files'] = False
-
+
h5_file.close()
@endcode
<h3>Resonance frequency and bandgap</h3>
-The results are analyzed in the
+The results are analyzed in the
[jupyter notebook](https://github.com/dealii/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb)
with the following code
@code{.py}
@f]
The square of the amplitude of the mechanical resonance $a^2$ as a function of the frequency
-has a gaussian shape
+has a gaussian shape
@f[
a^2 = a_\textrm{max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
@f]
<h3>Mode profile</h3>
-We can inspect the mode profile with Paraview or Visit.
+We can inspect the mode profile with Paraview or Visit.
As we have discussed, at resonance all the mechanical
energy is transmitted and the amplitude of motion is amplified inside the cavity.
It can be observed that the PMLs are quite effective to truncate the solution.
origin removed. In addition, we use $\varepsilon=0.005$, $\boldsymbol{\beta} =
[-\sin(\pi/6),\,\cos(\pi/6)]$, $f=0$, and Dirichlet boundary values
@f{align*}{
-g = \left\{\begin{array}{ll} 1 & \text{if } x=-1 \text{ or } y=-1,\,x\geq 0.5 \\
+g = \left\{\begin{array}{ll} 1 & \text{if } x=-1 \text{ or } y=-1,\,x\geq 0.5 \\
0 & \text{otherwise} \end{array}\right.
@f}
<th>48</th>
<td>3</th>
<td>3</th>
- <td>3</th>
+ <td>3</th>
<th></th>
<td>3</th>
<td>3</th>
- <td>3</th>
+ <td>3</th>
</tr>
<tr>
<th>128</th>
<td>6</th>
<td>6</th>
<td>6</th>
- <th></th>
+ <th></th>
+ <td>6</th>
<td>6</th>
<td>6</th>
- <td>6</th>
</tr>
<tr>
<th>512</th>
<th></th>
<td>9</th>
<td>9</th>
- <td>9</th>
+ <td>9</th>
</tr>
<tr>
<th>2048</th>
<th>2176</th>
<td>15</th>
<td>15</th>
- <td>15</th>
+ <td>15</th>
<th></th>
<td>13</th>
<td>13</th>
- <td>13</th>
+ <td>13</th>
</tr>
<tr>
<th>8192</th>
<th>8448</th>
<td>18</th>
<td>18</th>
- <td>18</th>
+ <td>18</th>
<th></th>
<td>15</th>
<td>15</th>
- <td>15</th>
+ <td>15</th>
</tr>
<tr>
<th>32768</th>
<th>33280</th>
<td>20</th>
<td>20</th>
- <td>20</th>
+ <td>20</th>
<th></th>
<td>16</th>
<td>16</th>
- <td>16</th>
+ <td>16</th>
</tr>
<tr>
<th>131072</th>
<th>132096</th>
<td>20</th>
<td>20</th>
- <td>20</th>
+ <td>20</th>
<th></th>
<td>16</th>
<td>16</th>
- <td>16</th>
+ <td>16</th>
</tr>
</table>
<th>48</th>
<td>2</th>
<td>2</th>
- <td>3</th>
+ <td>3</th>
<th></th>
<td>2</th>
<td>2</th>
- <td>3</th>
+ <td>3</th>
</tr>
<tr>
<th>128</th>
<td>5</th>
<td>5</th>
<td>7</th>
- <th></th>
+ <th></th>
<td>5</th>
<td>5</th>
- <td>7</th>
+ <td>7</th>
</tr>
<tr>
<th>512</th>
<th></th>
<td>7</th>
<td>7</th>
- <td>12</th>
+ <td>12</th>
</tr>
<tr>
<th>2048</th>
<td>10</th>
<td>12</th>
<td>15</th>
- <th></th>
+ <th></th>
<td>8</th>
<td>10</th>
- <td>17</th>
+ <td>17</th>
</tr>
<tr>
<th>8192</th>
<th>8448</th>
<td>11</th>
<td>15</th>
- <td>19</th>
+ <td>19</th>
<th></th>
<td>10</th>
<td>11</th>
- <td>20</th>
+ <td>20</th>
</tr>
<tr>
<th>32768</th>
<th>33280</th>
<td>12</th>
<td>16</th>
- <td>20</th>
+ <td>20</th>
<th></th>
<td>10</th>
<td>12</th>
- <td>21</th>
+ <td>21</th>
</tr>
<tr>
<th>131072</th>
<th>132096</th>
<td>12</th>
<td>16</th>
- <td>19</th>
+ <td>19</th>
<th></th>
<td>11</th>
<td>12</th>
- <td>21</th>
+ <td>21</th>
</tr>
</table>
iteration. Block smoothers here are at a disadvantage as they are
having to construct and invert a cell matrix for each cell. Here is a
comparison of solve times for a $Q_3$ element with 74,496 DoFs:
-
+
<table align="center" class="doxtable">
<tr>
<th colspan="1">$Q_3$</th>
@code
IndexSet locally_owned_dofs, locally_relevant_dofs;
...fill the two IndexSet objects...
-
+
// Create the ReadWriteVector using an IndexSet instead of the size
LinearAlgebra::ReadWriteVector<double> owned_rw_vector(locally_owned_dofs);
which implies that they overlap between different MPI
processes. Consequently, the elements stored in that vector on one
machine may not coincide with the ones stored by the GPU on that
-machine, requiring MPI communication to import them.
+machine, requiring MPI communication to import them.
In all of these cases, while importing a vector, values can either be
inserted (using VectorOperation::insert) or added to prior content of
optimal solver would have a number of iterations that is independent
of the size of the problem. But having such a solver would require
using a better preconditioner than the identity matrix we have used here.
-
+
<a name="extensions"></a>
<h3> Possible extensions </h3>
would not produce reasonable new points. In fact, defining a
manifold description for things that are supposed "to look round"
but might have points at or close to the origin is surprisingly very
- difficult.
+ difficult.
- A similar thing happens at the origin
of the three-dimensional ball when one tries to attach a spherical manifold to
the whole volume – in this case, the computation of new manifold points
is evaluated on a given cell
subject to this manifold. Using this manifold on the shaded cells of the
coarse grid of the disk (i.e., not only in the outer-most layer of
-cells) produces the following mesh upon four global
+cells) produces the following mesh upon four global
steps of refinement:
@image html circular_mesh_transfinite_interpolation.png ""
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
<li> Non-homogeneous Neumann boundary conditions for the Helmholtz equation.
</ol>
Besides these topics, again a variety of improvements and tricks will be
-shown.
+shown.
<h3>Verification of correctness</h3>
\\
{| u-u_h |}_{H^1(K)} &=& \left( \int_K |\nabla(u-u_h)|^2 \; dx \right)^{1/2},
\\
- {\| u-u_h \|}_{H^1(K)} &=& \left( {\| u-u_h \|}^2_{L_2(K)}
+ {\| u-u_h \|}_{H^1(K)} &=& \left( {\| u-u_h \|}^2_{L_2(K)}
+{| u-u_h |}^2_{H^1(K)} \right)^{1/2}.
@f}
(All these norms and semi-norms can also be evaluated with weighting functions,
rule with points only at points where super-convergence might occur, such as
the Gauss points of the lowest-order Gauss quadrature formula for which the
integrals in the assembly of the matrix is correct (e.g., for linear elements,
-do not use the QGauss(2) quadrature formula). In fact, this is generally good
+do not use the QGauss(2) quadrature formula). In fact, this is generally good
advice also for the other norms: if your quadrature points are fortuitously
chosen at locations where the error happens to be particularly small due to
superconvergence, the computed error will look like it is much smaller than
Because we want to verify the convergence of our numerical solution $u_h$,
we want a setup so that we know the exact solution $u$. This is where
-the Method of Manufactured Solutions comes in. To this end, let us
+the Method of Manufactured Solutions comes in. To this end, let us
choose a function
@f[
\bar u(x) = \sum_{i=1}^3 \exp\left(-\frac{|x-x_i|^2}{\sigma^2}\right)
@f]
-where the centers $x_i$ of the exponentials are
+where the centers $x_i$ of the exponentials are
$x_1=(-\frac 12,\frac 12)$,
$x_2=(-\frac 12,-\frac 12)$, and
$x_3=(\frac 12,-\frac 12)$,
Number of active cells: 15196
Number of degrees of freedom: 15912
-cycle cells dofs L2 H1 Linfty
- 0 64 81 1.840e+00 2.858e+00 1.835e+00
- 1 124 157 5.190e-02 1.200e+00 1.344e-01
- 2 280 341 1.439e-02 7.892e-01 7.554e-02
- 3 577 690 8.627e-03 5.061e-01 2.805e-02
- 4 1099 1264 3.217e-03 3.030e-01 1.073e-02
- 5 2191 2452 1.445e-03 2.097e-01 5.073e-03
- 6 4165 4510 8.387e-04 1.460e-01 2.013e-03
- 7 7915 8440 7.051e-04 1.053e-01 1.804e-03
- 8 15196 15912 2.774e-04 7.463e-02 6.911e-04
+cycle cells dofs L2 H1 Linfty
+ 0 64 81 1.840e+00 2.858e+00 1.835e+00
+ 1 124 157 5.190e-02 1.200e+00 1.344e-01
+ 2 280 341 1.439e-02 7.892e-01 7.554e-02
+ 3 577 690 8.627e-03 5.061e-01 2.805e-02
+ 4 1099 1264 3.217e-03 3.030e-01 1.073e-02
+ 5 2191 2452 1.445e-03 2.097e-01 5.073e-03
+ 6 4165 4510 8.387e-04 1.460e-01 2.013e-03
+ 7 7915 8440 7.051e-04 1.053e-01 1.804e-03
+ 8 15196 15912 2.774e-04 7.463e-02 6.911e-04
Solving with Q1 elements, global refinement
===========================================
Number of active cells: 16384
Number of degrees of freedom: 16641
-cycle cells dofs L2 H1 Linfty
- 0 64 81 1.840e+00 2.858e+00 1.835e+00
- 1 256 289 3.570e-02 1.199e+00 1.307e-01
- 2 1024 1089 1.192e-02 7.565e-01 7.168e-02
- 3 4096 4225 3.047e-03 3.823e-01 2.128e-02
- 4 16384 16641 7.660e-04 1.917e-01 5.554e-03
+cycle cells dofs L2 H1 Linfty
+ 0 64 81 1.840e+00 2.858e+00 1.835e+00
+ 1 256 289 3.570e-02 1.199e+00 1.307e-01
+ 2 1024 1089 1.192e-02 7.565e-01 7.168e-02
+ 3 4096 4225 3.047e-03 3.823e-01 2.128e-02
+ 4 16384 16641 7.660e-04 1.917e-01 5.554e-03
-n cells H1 L2
-0 64 2.858e+00 - - 1.840e+00 - -
-1 256 1.199e+00 2.38 1.25 3.570e-02 51.54 5.69
-2 1024 7.565e-01 1.58 0.66 1.192e-02 2.99 1.58
-3 4096 3.823e-01 1.98 0.98 3.047e-03 3.91 1.97
-4 16384 1.917e-01 1.99 1.00 7.660e-04 3.98 1.99
+n cells H1 L2
+0 64 2.858e+00 - - 1.840e+00 - -
+1 256 1.199e+00 2.38 1.25 3.570e-02 51.54 5.69
+2 1024 7.565e-01 1.58 0.66 1.192e-02 2.99 1.58
+3 4096 3.823e-01 1.98 0.98 3.047e-03 3.91 1.97
+4 16384 1.917e-01 1.99 1.00 7.660e-04 3.98 1.99
Solving with Q2 elements, global refinement
===========================================
Number of active cells: 16384
Number of degrees of freedom: 66049
-cycle cells dofs L2 H1 Linfty
- 0 64 289 1.606e-01 1.278e+00 3.029e-01
- 1 256 1089 7.638e-03 5.248e-01 4.816e-02
- 2 1024 4225 8.601e-04 1.086e-01 4.827e-03
- 3 4096 16641 1.107e-04 2.756e-02 7.802e-04
- 4 16384 66049 1.393e-05 6.915e-03 9.971e-05
+cycle cells dofs L2 H1 Linfty
+ 0 64 289 1.606e-01 1.278e+00 3.029e-01
+ 1 256 1089 7.638e-03 5.248e-01 4.816e-02
+ 2 1024 4225 8.601e-04 1.086e-01 4.827e-03
+ 3 4096 16641 1.107e-04 2.756e-02 7.802e-04
+ 4 16384 66049 1.393e-05 6.915e-03 9.971e-05
-n cells H1 L2
-0 64 1.278e+00 - - 1.606e-01 - -
-1 256 5.248e-01 2.43 1.28 7.638e-03 21.03 4.39
-2 1024 1.086e-01 4.83 2.27 8.601e-04 8.88 3.15
-3 4096 2.756e-02 3.94 1.98 1.107e-04 7.77 2.96
-4 16384 6.915e-03 3.99 1.99 1.393e-05 7.94 2.99
+n cells H1 L2
+0 64 1.278e+00 - - 1.606e-01 - -
+1 256 5.248e-01 2.43 1.28 7.638e-03 21.03 4.39
+2 1024 1.086e-01 4.83 2.27 8.601e-04 8.88 3.15
+3 4096 2.756e-02 3.94 1.98 1.107e-04 7.77 2.96
+4 16384 6.915e-03 3.99 1.99 1.393e-05 7.94 2.99
Solving with Q2 elements, adaptive refinement
===========================================
Number of active cells: 14212
Number of degrees of freedom: 64731
-cycle cells dofs L2 H1 Linfty
- 0 64 289 1.606e-01 1.278e+00 3.029e-01
- 1 124 577 7.891e-03 5.256e-01 4.852e-02
- 2 289 1353 1.070e-03 1.155e-01 4.868e-03
- 3 547 2531 5.962e-04 5.101e-02 1.876e-03
- 4 1057 4919 1.977e-04 3.094e-02 7.923e-04
- 5 2059 9223 7.738e-05 1.974e-02 7.270e-04
- 6 3913 17887 2.925e-05 8.772e-03 1.463e-04
- 7 7441 33807 1.024e-05 4.121e-03 8.567e-05
- 8 14212 64731 3.761e-06 2.108e-03 2.167e-05
+cycle cells dofs L2 H1 Linfty
+ 0 64 289 1.606e-01 1.278e+00 3.029e-01
+ 1 124 577 7.891e-03 5.256e-01 4.852e-02
+ 2 289 1353 1.070e-03 1.155e-01 4.868e-03
+ 3 547 2531 5.962e-04 5.101e-02 1.876e-03
+ 4 1057 4919 1.977e-04 3.094e-02 7.923e-04
+ 5 2059 9223 7.738e-05 1.974e-02 7.270e-04
+ 6 3913 17887 2.925e-05 8.772e-03 1.463e-04
+ 7 7441 33807 1.024e-05 4.121e-03 8.567e-05
+ 8 14212 64731 3.761e-06 2.108e-03 2.167e-05
@endcode
To see this, create a plot in log-log style with the number of unknowns on the
x axis and the L2 error on the y axis. You can add reference lines for
$h^2=N^{-1}$ and $h^3=N^{-3/2}$ and check that global and adaptive refinement
-follow those.
+follow those.
Note that changing the half width of the peaks influences if adaptive or
global refinement is more efficient (if the solution is very smooth, local
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
-\sum_{e \in \mathbb{F}} \int_{e}
\average{\frac{\partial^2 v_h}{\partial \mathbf n^2}}
\jump{\frac{\partial u_h}{\partial \mathbf n}}
-+ \sum_{e \in \mathbb{F}}
++ \sum_{e \in \mathbb{F}}
\frac{\sigma}{h_e}\int_e
\jump{\frac{\partial v_h}{\partial \mathbf n}}
\jump{\frac{\partial u_h}{\partial \mathbf n}}.
\average{\frac{\partial^2 v_h}{\partial \mathbf n^2}}
\jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds
\\
-&+ \sum_{e \in \mathbb{F}}
+&+ \sum_{e \in \mathbb{F}}
\frac{\gamma}{h_e}
\int_e
\jump{\frac{\partial v_h}{\partial \mathbf n}}
\jump{\frac{\partial u_h}{\partial \mathbf n}} \ ds,
@f}
-and
+and
@f{align*}{
\mathcal{F}(v_h)&:=\sum_{K \in \mathbb{T}}\int_{K} v_h f \ dx
+
this, but we also verify that this conjecture is experimentally true
below.
-
-
+
+
<h3>Other Boundary Conditions</h3>
We remark that the derivation of the $C^0$ IP method for the
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
-# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
+# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
* not write any data at all, the process should call the function
* DataSet::write_none() because the operation is *collective* and all the MPI
* processes have to contribute to the call, even if they don't have data to
- * write.
+ * write.
* @code
* std::vector<hsize_t> dataset_dimensions = {50, 30};
* auto dataset = group.create_dataset<double>("name", dataset_dimensions);
* CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> *fe_eval,
* const unsigned int q_point) const;
* \endcode
- *
+ *
* @deprecated Use apply_for_each_quad_point() instead.
*/
// clang-format on
// clang-format off
/**
- * Same as above, except that the functor @func only takes a single input
+ * Same as above, except that the functor @func only takes a single input
* argument (fe_eval) and computes the quadrature point from the thread id.
*
* @p func needs to define
projects / dealii.git / commitdiff