# Boilerplate: The Make.global_options expects variables to be set to
# yes, as is common for Makefiles.
#
- COND_SET_TO_YES(DEAL_II_WITH_TBB MAKEFILE_enablethreads)
+ COND_SET_TO_YES(DEAL_II_WITH_THREADS MAKEFILE_enablethreads)
COND_SET_TO_YES(DEAL_II_WITH_FUNCTIONPARSER MAKEFILE_enableparser)
COND_SET_TO_YES(BUILD_SHARED_LIBS MAKEFILE_enableshared)
OBJEXT = o
EXEEXT =
+shared-lib-suffix = .so
+static-lib-suffix = .a
+
+ifeq ($(enable-shared),yes)
+ lib-suffix = $(shared-lib-suffix)
+else
+ lib-suffix = $(static-lib-suffix)
+endif
+
# set paths to all the libraries we need:
lib-deal2.o = @MAKEFILE_TARGETS_RELEASE@
STRING(REGEX REPLACE
"^[0-9]+\\.[0-9]+\\.([0-9]+).*$" "\\1"
P4EST_VERSION_SUBMINOR "${P4EST_VERSION}")
+
+ # Now for the patch number such as in 0.3.4.1. If there
+ # is no patch number, then the REGEX REPLACE will fail,
+ # setting P4EST_VERSION_PATCH to P4EST_VERSION. If that
+ # is the case, then set the patch number to zero
STRING(REGEX REPLACE
- "^[0-9]+\\.[0-9]+\\.[0-9]+\\.([0-9]+).*$" "\\1"
+ "^[0-9]+\\.[0-9]+\\.[0-9]+\\.([0-9]+)?.*$" "\\1"
P4EST_VERSION_PATCH "${P4EST_VERSION}")
+ IF(${P4EST_VERSION_PATCH} STREQUAL "${P4EST_VERSION}")
+ SET(P4EST_VERSION_PATCH "0")
+ ENDIF()
+
MARK_AS_ADVANCED(P4EST_DIR)
ELSE()
2.5 format.
<li><em>Matthias Maier:</em>
- CMake build system. Periodic boundary conditions. Enhancements throughout the library.
+ CMake build system for the library and the testsuite. Periodic
+ boundary conditions. Enhancements throughout the library.
<li><em>Cataldo Manigrasso:</em>
Work on the codimension-one meshes, DoFHandler, and finite
Packaging and configuration issues.
<li><em>Bruno Turcksin:</em>
- Extending deal.ii for 64-bit integer support, various other changes.
+ Extending deal.ii for 64-bit integer support. Converting the
+ testsuite to CMake. Various other changes.
<li><em>Kainan Wang:</em>
Extending deal.ii for 64-bit integer support.
<code>-NOTFOUND</code> and may be set by hand.
</p>
+ <h4>Library conflicts</h4>
+ <p>Caveat: if you have a set of standard libraries in the default
+ location, say <code>/usr/lib</code> and a set of
+ private versions of the same libraries, for instance because you
+ need different revisions sometimes, in your own library directory,
+ you may receive an error message of the form:</p>
+
+<pre>
+CMake Warning at source/CMakeLists.txt:65 (ADD_LIBRARY):
+ Cannot generate a safe runtime search path for target deal_II.g because
+ files in some directories may conflict with libraries in implicit
+ directories:
+
+ runtime library [libtbb.so.2] in /usr/lib may be hidden by files in:
+ /my/private/lib
+
+ Some of these libraries may not be found correctly.
+</pre>
+
+<p>This is not a problem of CMake or deal.II, but rather a general
+ Linux problem. In order to fix this, you have two options:</p>
+<ol>
+<li>Choose all libraries either from your private directory or from
+ the standard one.</li>
+<li>Install all your private library versions in different directories.</li>
+</ol>
+
<a name="configureoverride"></a>
<h3> Manual override </h3>
<p><b>Warning:</b> Do not do this unless absolutely necessary!</p>
<h2>The regression tests</h2>
<p>
- deal.II has a testsuite that, at the time this article is written (mid-2011),
- has some 2300 small programs (growing by roughly one per day) that we run
+ deal.II has a testsuite that, at the time this article is written (mid-2013),
+ has some 2,900 small programs (growing by roughly one per day) that we run
every time we make a
change to make sure that no existing functionality is broken. The
expected output is also stored in our subversion archive, and when you
in both cases to make sure that future changes do not break what
we have just checked in. In addition, some machines run the tests
every night and send the results back home; this is then converted
- into <a href="http://www.dealii.org/cgi-bin/regression_quick.pl"
+ into <a href="http://dealii.mathsim.eu/cgi-bin/regression_quick.pl"
target="body">a webpage showing the status of our regression
tests</a>.
</p>
line like
<pre>
- cmake -DP4EST_DIR=/path/to/installation/FAST -DDEAL_II_WITH_P4EST=ON -DDEAL_II_WITH_MPI=ON <...>
+ cmake -DP4EST_DIR=/path/to/installation -DDEAL_II_WITH_P4EST=ON -DDEAL_II_WITH_MPI=ON <...>
</pre>
if the p4est library isn't picked up automatically. Note the presence
of <code>/FAST</code> at the end of the path necessary when using
<ol>
<li>
- Fixed: The ParameterHandler class can now deal with including one parameter
+ New: All vector classes now have functions <code>extract_subvector_to()</code>
+ that allow extracting not just a single value but a whole set.
+ <br>
+ (Fahad Alrasched, 2013/09/02)
+ </li>
+
+ <li>
+ Fixed: <code>common/Make.global_options</code> now exports enable-threads
+ correctly, furthermore, <code>lib-suffix</code>, <code>shared-lib-suffix</code>
+ and <code>static-lib-suffix</code> are now exported as well for better legacy
+ support.
+ <br>
+ (Matthias Maier, 2013/08/30)
+ </li>
+
+ <li>
+ New: The ParameterHandler class can now deal with including one parameter
file from another.
<br>
(Wolfgang Bangerth, 2013/08/25)
// will generate Gnuplot output, which consists of a set of lines
// describing the mapped triangulation. By default, only one line
// is drawn for each face of the triangulation, but since we want
- // to explicitely see the effect of the mapping, we want to have
+ // to explicitly see the effect of the mapping, we want to have
// the faces in more detail. This can be done by passing the
// output object a structure which contains some flags. In the
// present case, since Gnuplot can only draw straight lines, we
solutions of Laplace's equation are only in $H^1$, which does not allow for
the definition of point values because it is not a subset of the continuous
functions. Therefore, even though fixing one node is allowed for
- discretitized functions, it is not for continuous functions, and one can
+ discretized functions, it is not for continuous functions, and one can
often see this in a resulting error spike at this point in the numerical
solution.
consider the objects describing mappings as a black box which you need not
worry about, because their only uses seem to be to be passed to places deep
inside the library where functions know how to handle them (i.e. in the
-<code>FEValues</code> classes and their descendents).
+<code>FEValues</code> classes and their descendants).
The tricky point in this program is the use of the mean value
constraint. Fortunately, there is a class in the library which knows how to
case to use a higher order mapping, not because it improves the order
of convergence but just to reduce the constant before the convergence
order. On the other hand, using a cubic mapping only improves the
-result further insignicantly, except for the case of very coarse
+result further insignificantly, except for the case of very coarse
grids.
// the mean value of the degrees of freedom on the boundary shall be
// zero. For this, we first want a list of those nodes which are actually
// at the boundary. The <code>DoFTools</code> class has a function that
- // returns an array of boolean values where <code>true</code> indicates
+ // returns an array of Boolean values where <code>true</code> indicates
// that the node is at the boundary. The second argument denotes a mask
// selecting which components of vector valued finite elements we want to
// be considered. This sort of information is encoded using the
In previous examples, you have seen how the library itself is broken
-up into several complexes each building atop the underying ones, but
+up into several complexes each building atop the underlying ones, but
relatively independent of the other ones:
<ol>
<li>the triangulation class complex, with associated iterator classes;
documentation, which would bring that to maybe 1200 lines. Yet, compared to
the applications listed above, this is still small, as they are 20 to
25 times as large. For programs as large, a proper design right from
-the start is thus indispensible. Otherwise, it will have to be
+the start is thus indispensable. Otherwise, it will have to be
redesigned at one point in its life, once it becomes too large to be
manageable.
and when to reload it again; if this is not done in an organized
fashion, then you end up with data released too early, loaded too
late, or not released at all). Although the present example program
-thus draws from sevelar years of experience, it is certainly not
+thus draws from several years of experience, it is certainly not
without flaws in its design, and in particular might not be suited for
an application where the objective is different. It should serve as an
inspiration for writing your own application in a modular way, to
// The second of this pair of functions takes a range of cell iterators,
// and assembles the system matrix on this part of the domain. Since it's
// actions have all been explained in previous programs, we do not comment
- // on it any more, except for one pointe below.
+ // on it any more, except for one point below.
template <int dim>
void
Solver<dim>::assemble_matrix (LinearSystem &linear_system,
// will always be unlocked when we exit this part of the program,
// whether the operation completed successfully or not, whether the
// exit path was something we implemented willfully or whether the
- // function was exited by an exception that we did not forsee.
+ // function was exited by an exception that we did not foresee.
//
// deal.II implements the scoped locking pattern in the
// Treads::Mutex::ScopedLock class: it takes the mutex in the
// @sect4{Local refinement by the Kelly error indicator}
// The second class implementing refinement strategies uses the Kelly
- // refinemet indicator used in various example programs before. Since this
+ // refinement indicator used in various example programs before. Since this
// indicator is already implemented in a class of its own inside the
// deal.II library, there is not much t do here except cal the function
// computing the indicator, then using it to select a number of cells for
</td>
</tr>
</table>
-Note the assymetry of the grids compared with those we obtained for
+Note the asymmetry of the grids compared with those we obtained for
the point evaluation, which is due to the directionality of the
x-derivative for which we tailored the refinement criterion.
// afterwards...
double point_derivative = 0;
- // ...then have some objects of which the meaning wil become clear
+ // ...then have some objects of which the meaning will become clear
// below...
QTrapez<dim> vertex_quadrature;
FEValues<dim> fe_values (dof_handler.get_fe(),
// We will then implement two such classes, for the evaluation of a point
// value and of the derivative of the solution at that point. For these
// functionals we already have the corresponding evaluation objects, so they
- // are comlementary.
+ // are complementary.
namespace DualFunctional
{
// @sect4{The DualFunctionalBase class}
// Next we have the function that is called to estimate the error on a
// subset of cells. The function may be called multiply if the library was
- // configured to use multi-threading. Here it goes:
+ // configured to use multithreading. Here it goes:
template <int dim>
void
WeightedResidual<dim>::
// Then calculate the start cell for this thread. We let the different
// threads run on interleaved cells, i.e. for example if we have 4
- // threads, then the first thread treates cells 0, 4, 8, etc, while the
+ // threads, then the first thread treats cells 0, 4, 8, etc, while the
// second threads works on cells 1, 5, 9, and so on. The reason is that
// it takes vastly more time to work on cells with hanging nodes than on
// regular cells, but such cells are not evenly distributed across the
// element solution at the quadrature points on the other side of the
// face, i.e. from the neighboring cell.
//
- // For this, do a sanity check before: make sure that the neigbor
+ // For this, do a sanity check before: make sure that the neighbor
// actually exists (yes, we should not have come here if the neighbor
// did not exist, but in complicated software there are bugs, so better
// check this), and if this is not the case throw an error.
// the name <code>neighbor_neighbor</code>:
const unsigned int
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
- // Then define an abbreviation for the neigbor cell, initialize the
+ // Then define an abbreviation for the neighbor cell, initialize the
// <code>FEFaceValues</code> object on that cell, and extract the
// gradients on that cell:
const active_cell_iterator neighbor = cell->neighbor(face_no);
// Then find out which neighbor the present cell is of the adjacent
// cell. Note that we will operate on the children of this adjacent
// cell, but that their orientation is the same as that of their mother,
- // i.e. the neigbor direction is the same.
+ // i.e. the neighbor direction is the same.
const unsigned int
neighbor_neighbor = cell->neighbor_of_neighbor (face_no);
right boundary values, the Newton update $\delta u^{n}$ should have zero
boundary conditions, in order to have the right boundary condition after
adding both. In the first Newton step, we are starting with the solution
-$u^{0}\equiv 0$, the Newton update still has to deliever the right boundary
+$u^{0}\equiv 0$, the Newton update still has to deliver the right boundary
condition to the solution $u^{1}$.
@f]
Using the basis functions as test functions and defining $a_{n}:=\frac{1}
-{\sqrt{1+|\nabla u^{n}|^{2}}}$, we can rewrite the weak formualtion:
+{\sqrt{1+|\nabla u^{n}|^{2}}}$, we can rewrite the weak formulation:
@f[
\sum_{j=0}^{N-1}\left[ \left( \nabla \varphi_{i} , a_{n} \nabla \varphi_{j} \right) -
\nabla u^{n}\right) \qquad \forall i=0,\dots ,N-1
@f]
-where the solution $\delta u^{n}$ is given by the coefficents $\delta U^{n}_{j}$.
+where the solution $\delta u^{n}$ is given by the coefficients $\delta U^{n}_{j}$.
This linear system of equations can be rewritten as:
@f[
It is clearly visible, that the solution minimizes the surface
after each refinement. The solution converges to a picture one
-would imagine a soapbubble to be that is located inside a wire loop
+would imagine a soap bubble to be that is located inside a wire loop
that is bent like
the boundary. Also it is visible, how the boundary
is smoothed out after each refinement. On the coarse mesh,
#include <fstream>
#include <iostream>
-// We will use adaptive mesh refinement between Newton interations. To do so,
+// We will use adaptive mesh refinement between Newton iterations. To do so,
// we need to be able to work with a solution on the new mesh, although it was
// computed on the old one. The SolutionTransfer class transfers the solution
// from the old to the new mesh:
#include <deal.II/numerics/solution_transfer.h>
-// We then open a namepsace for this program and import everything from the
+// We then open a namespace for this program and import everything from the
// dealii namespace into it, as in previous programs:
namespace Step15
{
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
-// These, now, are the include necessary for the multi-level methods. The
+// These, now, are the include necessary for the multilevel methods. The
// first two declare classes that allow us to enumerate degrees of freedom not
// only on the finest mesh level, but also on intermediate levels (that's what
// the MGDoFHandler class does) as well as allow to access this information
// Now for the things that concern the multigrid data structures. First,
- // we resize the multi-level objects to hold matrices and sparsity
+ // we resize the multilevel objects to hold matrices and sparsity
// patterns for every level. The coarse level is zero (this is mandatory
// right now but may change in a future revision). Note that these
// functions take a complete, inclusive range here (not a starting index
// integration core is the same as above, but the loop below will go over
// all existing cells instead of just the active ones, and the results must
// be entered into the correct matrix. Note also that since we only do
- // multi-level preconditioning, no right-hand side needs to be assembled
+ // multilevel preconditioning, no right-hand side needs to be assembled
// here.
//
// Before we go there, however, we have to take care of a significant amount
// that are located on interfaces between adaptively refined levels, and
// those that lie on the interface but also on the exterior boundary of
// the domain. As in many other parts of the library, we do this by using
- // boolean masks, i.e. vectors of booleans each element of which indicates
+ // Boolean masks, i.e. vectors of Booleans each element of which indicates
// whether the corresponding degree of freedom index is an interface DoF
// or not. The <code>MGConstraints</code> already computed the information
// for us when we called initialize in <code>setup_system()</code>.
MGTransferPrebuilt<Vector<double> > mg_transfer(hanging_node_constraints, mg_constrained_dofs);
// Now the prolongation matrix has to be built. This matrix needs to take
// the boundary values on each level into account and needs to know about
- // the indices at the refinement egdes. The <code>MGConstraints</code>
+ // the indices at the refinement edges. The <code>MGConstraints</code>
// knows about that so pass it as an argument.
mg_transfer.build_matrices(mg_dof_handler);
// to the present process, but then we need to distribute the refinement
// indicators into a distributed vector so that all processes have the
// values of the refinement indicator for all cells. But then, in order for
- // each process to refine its copy of the mesh, they need to have acces to
+ // each process to refine its copy of the mesh, they need to have access to
// all refinement indicators locally, so they have to copy the global vector
// back into a local one. That's a little convoluted, but thinking about it
// quite straightforward nevertheless. So here's how we do it:
// <code>local_error_per_cell</code> vector. The elements of this vector
// for cells not on the present process are zero. However, since all
// processes have a copy of a copy of the entire triangulation and need to
- // keep these copies in synch, they need the values of refinement
+ // keep these copies in sync, they need the values of refinement
// indicators for all cells of the triangulation. Thus, we need to
// distribute our results. We do this by creating a distributed vector
// where each process has its share, and sets the elements it has
corresponding stress update. However, the rotational component does not
necessarily induce a nonzero stress update (think, in 2d, for example of the
situation where $\Delta\mathbf{u}=(y, -x)^T$, with which $\varepsilon(\Delta
-\mathbf{u})=0$). Nevertheless, if the the material was pre-stressed in a certain
+\mathbf{u})=0$). Nevertheless, if the the material was prestressed in a certain
direction, then this direction will be rotated along with the material. To
this end, we have to define a rotation matrix $R(\Delta \mathbf{u}^n)$ that
describes, in each point the rotation due to the displacement increments. It
// the system, direct what has to be solved in each time step, a function
// that solves the linear system that arises in each timestep (and returns
// the number of iterations it took), and finally output the solution
- // vector on the currect mesh:
+ // vector on the correct mesh:
void create_coarse_grid ();
void setup_system ();
// @sect4{The public interface}
- // The next step is the definition of constructors and descructors. There
+ // The next step is the definition of constructors and destructors. There
// are no surprises here: we choose linear and continuous finite elements
// for each of the <code>dim</code> vector components of the solution, and a
// Gaussian quadrature formula with 2 points in each coordinate
// The next function is the one that sets up the data structures for a given
// mesh. This is done in most the same way as in step-17: distribute the
// degrees of freedom, then sort these degrees of freedom in such a way that
- // each processor gets a contiguous chunk of them. Note that subdivions into
+ // each processor gets a contiguous chunk of them. Note that subdivisions into
// chunks for each processor is handled in the functions that create or
// refine grids, unlike in the previous example program (the point where
// this happens is mostly a matter of taste; here, we chose to do it when
// adaptively refined).
//
// With this data structure, we can then go to the PETSc sparse matrix and
- // tell it to pre-allocate all the entries we will later want to write to:
+ // tell it to preallocate all the entries we will later want to write to:
system_matrix.reinit (mpi_communicator,
sparsity_pattern,
local_dofs_per_process,
// <code>cell-@>vertex_dof_index(v,d)</code> function that returns the index
// of the <code>d</code>th degree of freedom at vertex <code>v</code> of the
// given cell. In the present case, displacement in the k-th coordinate
- // direction corresonds to the kth component of the finite element. Using a
+ // direction corresponds to the k-th component of the finite element. Using a
// function like this bears a certain risk, because it uses knowledge of the
// order of elements that we have taken together for this program in the
// <code>FESystem</code> element. If we decided to add an additional
// variable, for example a pressure variable for stabilization, and happened
// to insert it as the first variable of the element, then the computation
- // below will start to produce non-sensical results. In addition, this
+ // below will start to produce nonsensical results. In addition, this
// computation rests on other assumptions: first, that the element we use
// has, indeed, degrees of freedom that are associated with vertices. This
// is indeed the case for the present Q1 element, as would be for all Qp
// To put this into larger perspective, we note that if we had previously
// available stresses in our model (which we assume do not exist for the
// purpose of this program), then we would need to interpolate the field of
- // pre-existing stresses to the quadrature points. Likewise, if we were to
+ // preexisting stresses to the quadrature points. Likewise, if we were to
// simulate elasto-plastic materials with hardening/softening, then we would
// have to store additional history variables like the present yield stress
// of the accumulated plastic strains in each quadrature
// displacement update so that the material in its new configuration
// accommodates for the difference between the external body and boundary
// forces applied during this time step minus the forces exerted through
- // pre-existing internal stresses. In order to have the pre-existing
+ // preexisting internal stresses. In order to have the preexisting
// stresses available at the next time step, we therefore have to update the
- // pre-existing stresses with the stresses due to the incremental
+ // preexisting stresses with the stresses due to the incremental
// displacement computed during the present time step. Ideally, the
// resulting sum of internal stresses would exactly counter all external
// forces. Indeed, a simple experiment can make sure that this is so: if we
// three matrices should be symmetric, it is not due to floating
// point round off: we get an asymmetry on the order of 1e-16 of
// the off-diagonal elements of the result. When assigning the
- // result to a <code>SymmetricTensor</code>, the constuctor of
+ // result to a <code>SymmetricTensor</code>, the constructor of
// that class checks the symmetry and realizes that it isn't
// exactly symmetric; it will then raise an exception. To avoid
// that, we explicitly symmetrize the result to make it exactly
subsection DX output parameters
- # A boolean field indicating whether neighborship information between cells
+ # A Boolean field indicating whether neighborship information between cells
# is to be written to the OpenDX output file
set Write neighbors = true
end
end
subsection Povray output parameters
- # Whether camera and lightling information should be put into an external
+ # Whether camera and lighting information should be put into an external
# file "data.inc" or into the POVRAY input file
set Include external file = true
// this program, I realized that there aren't all that many parameters this
// program can usefully ask for, or better, it turned out: declaring and
// querying these parameters was already done centralized in one place of
- // the libray, namely the <code>DataOutInterface</code> class that handles
+ // the library, namely the <code>DataOutInterface</code> class that handles
// exactly this -- managing parameters for input and output.
//
// So the second function call in this function is to let the
// with the subsection name, I like to use curly braces to force my editor
// to indent everything that goes into this sub-section by one level of
// indentation. In this sub-section, we shall have two entries, one that
- // takes a boolean parameter and one that takes a selection list of
+ // takes a Boolean parameter and one that takes a selection list of
// values, separated by the '|' character:
prm.enter_subsection ("Dummy subsection");
{
// Finally, let us note that if we were interested in the values
// of the parameters declared above in the dummy subsection, we
// would write something like this to extract the value of the
- // boolean flag (the <code>prm.get</code> function returns the
+ // Boolean flag (the <code>prm.get</code> function returns the
// value of a parameter as a string, whereas the
// <code>prm.get_X</code> functions return a value already
// converted to a different type):
to indicate two slightly different, but related things. The first is that we'd
like to represent the finite element solution as a linear combination of shape
function, in the form $u_h(\mathbf x) = \sum_{j=0}^{N-1} U_j \varphi_j(\mathbf
-x)$. Here, $U_j$ is a vector of expension coefficients. Because we don't know
+x)$. Here, $U_j$ is a vector of expansion coefficients. Because we don't know
their values yet (we will compute them as the solution of a linear or
nonlinear system), they are called "unknowns" or "degrees of freedom". The
second meaning of the term con be explained as follows: A mathematical
number of nonzero entries per row is five, and therefore independent of the
total size of the matrix.) Sparsity is one of the distinguishing feature of
the finite element method compared to, say, approximating the solution of a
-partial differential equation using a Taylor expension and matching
+partial differential equation using a Taylor expansion and matching
coefficients, or using a Fourier basis.
In practical terms, it is the sparsity of matrices that enables us to solve
patterns. We can visualize them using <a
href="http://www.gnuplot.info/">GNUPLOT</a> (one of the simpler visualization
programs; maybe not the easiest to use since it is command line driven, but
-also universally available on all linux and other unix-like systems):
+also universally available on all Linux and other Unix-like systems):
@code
examples/\step-2> gnuplot
// least as long as that of the <code>DoFHandler</code>; one way to make sure
// this is so is to make it static as well, in order to prevent its preemptive
// destruction. (However, the library would warn us if we forgot about this
-// and abort the program if that occured. You can check this, if you want, by
+// and abort the program if that occurred. You can check this, if you want, by
// removing the 'static' declaration.)
void distribute_dofs (DoFHandler<2> &dof_handler)
{
the x-velocity has to increase to get the flow through the narrow part,
something that can easily be seen in the left image. The middle image
represents inward flow in y-direction at the left end of the domain, and
-outward flow in y-directino at the right end of the domain.
+outward flow in y-direction at the right end of the domain.
Realistic flow computations for ground water or oil reservoir simulations will
not use a constant permeability. Here's a first, rather simple way to change
this situation: we use a permeability that decays very rapidly away from a
-central flowline until it hits a background value of 0.001. This is to mimick
+central flowline until it hits a background value of 0.001. This is to mimic
the behavior of fluids in sandstone: in most of the domain, the sandstone is
-homogenous and, while permeably to fluids, not overly so; on the other stone,
+homogeneous and, while permeable to fluids, not overly so; on the other stone,
the stone has cracked, or faulted, along one line, and the fluids flow much
-easier along this large crask. Here is how we could implement something like
+easier along this large crack. Here is how we could implement something like
this:
@code
template <int dim>
// binding simpler elements together into one larger element. In the present
// case, we want to couple a single RT(degree) element with a single
// DQ(degree) element. The constructor to <code>FESystem</code> that does
- // this requires us to specity first the first base element (the
+ // this requires us to specify first the first base element (the
// <code>FE_RaviartThomas</code> object of given degree) and then the number
// of copies for this base element, and then similarly the kind and number
- // of <code>FE_DGQ</code> elements. Note that the Raviart Thomas element
+ // of <code>FE_DGQ</code> elements. Note that the Raviart-Thomas element
// already has <code>dim</code> vector components, so that the coupled
// element will have <code>dim+1</code> vector components, the first
// <code>dim</code> of which correspond to the velocity variable whereas the
// used here, the Gauss points happen to be superconvergence points in
// which the pointwise error is much smaller (and converges with higher
// order) than anywhere else. These are therefore not particularly good
- // points for ingration. To avoid this problem, we simply use a
+ // points for integration. To avoid this problem, we simply use a
// trapezoidal rule and iterate it <code>degree+2</code> times in each
// coordinate direction (again as explained in step-7):
QTrapez<1> q_trapez;
// The main function we stole from step-6 instead of step-4. It is almost
// equal to the one in step-6 (apart from the changed class names, of course),
// the only exception is that we pass the degree of the finite element space
-// to the constructor of the mixed laplace problem (here, we use zero-th order
+// to the constructor of the mixed Laplace problem (here, we use zero-th order
// elements).
int main ()
{
Here, $p=p(\mathbf x, t), S=S(\mathbf x, t)$ are now time dependent
functions: while at every time instant the flow field is in
equilibrium with the pressure (i.e. we neglect dynamic
-accellerations), the saturation is transported along with the flow and
+accelerations), the saturation is transported along with the flow and
therefore changes over time, in turn affected the flow field again
through the dependence of the first equation on $S$.
because they consist of a network of interconnected pores. However, this
microscopic scale is out of the range of today's computer simulations, and we
have to be content with simulating things on the scale of meters. On that
-scale, however, fluid transport typically happens through a network of crasks
-in the rock, rather than throuh pores. However, cracks often result from
+scale, however, fluid transport typically happens through a network of cracks
+in the rock, rather than through pores. However, cracks often result from
external stress fields in the rock layer (for example from tectonic faulting)
and the cracks are therefore roughly aligned. This leads to a situation where
the permeability is often orders of magnitude larger in the direction parallel
"step-20":
<ul>
<li>A function that models a single, winding crack that snakes through the
- domain. In analgy to step-20, but taking care of the slightly
+ domain. In analogy to step-20, but taking care of the slightly
different geometry we have here, we describe this by the following function:
@f[
k(\mathbf x)
</ol>
Despite this wealth of possibilities, we will only use Dirichlet and
-(homogenous) Neumann boundary conditions in this tutorial program.
+(homogeneous) Neumann boundary conditions in this tutorial program.
<h3>Discretization</h3>
represent an exact inverse of the matrix $A$. Consequently, preconditioning
with the ILU will still require more than one iteration, unlike
preconditioning with the sparse direct solver. The inner solver will therefore
-take more time when multiplying with the Schur complement, a tradeoff
+take more time when multiplying with the Schur complement, a trade-off
unavoidable.
</ol>
</a>
Frequently, a sparse matrix contains a substantial amount of elements that
-acutally are zero when we are about to start a linear solve. Such elements are
+actually are zero when we are about to start a linear solve. Such elements are
introduced when we eliminate constraints or implement Dirichlet conditions,
where we usually delete all entries in constrained rows and columns, i.e., we
set them to zero. The fraction of elements that are present in the sparsity
consumption. Likewise, the implementation of the decomposition step in the
SparseILU class was very inefficient and has been replaced by one that is
about 10 times faster. Even the vmult function of the SparseILU has been
-improved to save about twenty precent of time. Small improvements were
+improved to save about twenty percent of time. Small improvements were
applied here and there. Moreover, the ConstraintMatrix object has been used
to eliminate a lot of entries in the sparse matrix that are eventually going
to be zero, see <a href="#constraint-matrix">the section on using advanced
steps in the program. Observe how the grid is refined in regions
where the solution rapidly changes: On the upper boundary, we have
Dirichlet boundary conditions that are -1 in the left half of the line
-and 1 in the right one, so there is an aprupt change at $x=0$. Likewise,
+and 1 in the right one, so there is an abrupt change at $x=0$. Likewise,
there are changes from Dirichlet to Neumann data in the two upper
corners, so there is need for refinement there as well:
<img src="http://www.dealii.org/images/steps/developer/step-22.3d.solution.png" alt="">
-The isocountours shown here as well are those of the pressure
+The isocontours shown here as well are those of the pressure
variable, showing the singularity at the point of discontinuous
velocity boundary conditions.
always around five, both in two and three dimensions.) To summarize, most work
is spent on solving linear systems with the same matrix $A$ over and over again.
What makes this look even worse is the fact that we
-actually invert a matrix that is about 95 precent the size of the total system
-matrix and stands for 85 precent of the non-zero entries in the sparsity
+actually invert a matrix that is about 95 percent the size of the total system
+matrix and stands for 85 percent of the non-zero entries in the sparsity
pattern. Hence, the natural question is whether it is reasonable to solve a
linear system with matrix $A$ for about 15 times when calculating the solution
to the block system.
Another idea to improve the situation even more would be to choose a
preconditioner that makes CG for the (0,0) matrix $A$ converge in a
mesh-independent number of iterations, say 10 to 30. We have seen such a
-canditate in step-16: multigrid.
+candidate in step-16: multigrid.
<h5>Block Schur complement preconditioner</h5>
Even with a good preconditioner for $A$, we still
vice versa.
<a href="http://en.wikipedia.org/wiki/GMRES#Comparison_with_other_solvers">Wikipedia</a>'s
article on the GMRES method gives a comparative presentation.
-A more comprehensive and well-founded comparsion can be read e.g. in the book by
+A more comprehensive and well-founded comparison can be read e.g. in the book by
J.W. Demmel (Applied Numerical Linear Algebra, SIAM, 1997, section 6.6.6).
For our specific problem with the ILU preconditioner for $A$, we certainly need
complicated than the examples shown above, but it is not hard to accommodate
for that.
-For example, by using the folllowing modification of the boundary values
+For example, by using the following modification of the boundary values
function
@code
template <int dim>
// As in step-20 and most other example programs, the next task is to define
// the data for the PDE: For the Stokes problem, we are going to use natural
- // boundary values on parts of the boundary (i.e. homogenous Neumann-type)
+ // boundary values on parts of the boundary (i.e. homogeneous Neumann-type)
// for which we won't have to do anything special (the homogeneity implies
// that the corresponding terms in the weak form are simply zero), and
// boundary conditions on the velocity (Dirichlet-type) on the rest of the
// Note that we initialize the triangulation with a MeshSmoothing argument,
// which ensures that the refinement of cells is done in a way that the
// approximation of the PDE solution remains well-behaved (problems arise if
- // grids are too unstructered), see the documentation of
+ // grids are too unstructured), see the documentation of
// <code>Triangulation::MeshSmoothing</code> for details.
template <int dim>
StokesProblem<dim>::StokesProblem (const unsigned int degree)
// require allocating more memory than necessary but isn't suitable for
// use as a basis of SparseMatrix or BlockSparseMatrix objects; in a
// second step we then copy this object into an object of
- // BlockSparsityPattern. This is entirely analgous to what we already did
+ // BlockSparsityPattern. This is entirely analogous to what we already did
// in step-11 and step-18.
//
// There is one snag again here, though: it turns out that using the
// After the discussion in the introduction and the definition of the
// respective classes above, the implementation of the <code>solve</code>
- // function is rather straigt-forward and done in a similar way as in
+ // function is rather straight-forward and done in a similar way as in
// step-20. To start with, we need an object of the
// <code>InverseMatrix</code> class that represents the inverse of the
// matrix A. As described in the introduction, the inverse is generated with
// An alternative that is cheaper to build, but needs more iterations
// afterwards, would be to choose a SSOR preconditioner with factor
// 1.2. It needs about twice the number of iterations, but the costs for
- // its generation are almost neglible.
+ // its generation are almost negligible.
SparseILU<double> preconditioner;
preconditioner.initialize (system_matrix.block(1,1),
SparseILU<double>::AdditionalData());
// objects we already know - so we simply multiply $p$ by $B^T$, subtract
// the right hand side and multiply by the inverse of $A$. At the end, we
// need to distribute the constraints from hanging nodes in order to
- // obtain a constistent flow field:
+ // obtain a consistent flow field:
{
system_matrix.block(0,1).vmult (tmp, solution.block(1));
tmp *= -1;
Both of these methods have advantages and disadvantages.
Traditionally, people have preferred the method of lines, since it
allows to use the very well developed machinery of high-order ODE
-solvers avaiable for the rather stiff ODEs resulting from this
+solvers available for the rather stiff ODEs resulting from this
approach, including step length control and estimation of the temporal
error.
is what we set the time step to. The fact that we set the time step and mesh
size individually in two different places is error prone: it is too easy to
refine the mesh once more but forget to also adjust the time step. @ref
-step_24 "step-24" shows a better way how to keep these things in synch.
+step_24 "step-24" shows a better way how to keep these things in sync.
<h3>The test case</h3>
sparsity_pattern.compress();
// Then comes a block where we have to initialize the 3 matrices we need
- // in the course of the program: the mass matrix, the laplace matrix, and
+ // in the course of the program: the mass matrix, the Laplace matrix, and
// the matrix $M+k^2\theta^2A$ used when solving for $U^n$ in each time
// step.
//
homogeneous acoustic properties. The basic acoustic generation equation in an
acoustically homogeneous medium can be described as follows: if $u$ is the
vector-valued displacement, then tissue certainly reacts to changes in
-pressure by accelleration:
+pressure by acceleration:
@f[
\rho \frac{\partial^2}{\partial t^2}u(t,\mathbf r) =
-\nabla p(t,\mathbf r).
results from the use of absorbing boundary conditions.
Above two equations can be rewritten in a matrix form with the pressure and its derivative as
-an unknown vecotor:
+an unknown vector:
@f[
\left(\begin{array}{cc}
M & -k\theta M \\
values evaluated at each detector location to disk. We then
draw them in plots. Experimental data were also collected for comparison.
Currently our experiments have only been done in two dimensions by
-circually scanning a single detector. The tissue sample here is a thin slice
+circularly scanning a single detector. The tissue sample here is a thin slice
in the X-Y plane (Z=0), and we assume that signals from other Z directions
won't contribute to the data. Consequently, we only have to compare
our experimental data with two dimensional simulated data.
circular mesh we use (see for example step-6 for a view of the
coarse mesh) is not isotropic either. The net result is that the signal fronts
are not sinusoidal unless the mesh is sufficiently fine. The right image is a
-lot better in this respect, though artificts in the form of trailing spurious
+lot better in this respect, though artifacts in the form of trailing spurious
waves can still be seen.
\f}
The iteration can be initialized with the old time step, i.e. $u^n_0 = u^{n-1}$,
and eventually it will produce a solution to the first equation of
-the split formulation (see above). For the time discretizaion of the
+the split formulation (see above). For the time discretization of the
sine-Gordon equation under consideration here, we have that
\f{eqnarray*}
F(u^n_l) &=& \left[ 1-k^2\theta^2\Delta \right] u^n_l -
letter the vector of coefficients (in the nodal basis) of a function
denoted by the same letter in lower case; e.g., $u^n = \sum_{i=1}^N
U^n_i \varphi_i$ where $U^n \in {R}^N$ and $u^n \in
-H^1(\Omega)$. Thus, the finite-dimensional version of the variational fomulation requires that we solve the following matrix equations at each time step:
+H^1(\Omega)$. Thus, the finite-dimensional version of the variational formulation requires that we solve the following matrix equations at each time step:
@f{eqnarray*}
F_h'(U^{n,l})\delta U^{n,l} &=& -F_h(U^{n,l}), \qquad
U^{n,l+1} = U^{n,l} + \delta U^{n,l}, \qquad U^{n,0} = U^{n-1}; \\
- k^2\theta^2N(u^n_l,u^{n-1})
\f}
Again, note that the first matrix equation above is, in fact, the
-defition of an iterative procedure, so it is solved multiple times
+definition of an iterative procedure, so it is solved multiple times
until a stopping criterion is met. Moreover, $M$ is the mass matrix,
i.e. $M_{ij} = \left( \varphi_i,\varphi_j \right)_{\Omega}$, $A$ is
the Laplace matrix, i.e. $A_{ij} = \left( \nabla \varphi_i, \nabla
// integrate these terms exactly. It is usually sufficient to just make sure
// that the right hand side is integrated up to the same order of accuracy
// as the discretization scheme is, but it may be possible to improve on the
- // constant in the asympotitic statement of convergence by choosing a more
+ // constant in the asymptotic statement of convergence by choosing a more
// accurate quadrature formula.
template <int dim>
void SineGordonProblem<dim>::compute_nl_term (const Vector<double> &old_data,
// @sect4{SineGordonProblem::compute_nl_matrix}
// This is the second function dealing with the nonlinear scheme. It
- // computes the matrix $N(\cdot,\cdot)$, whicih appears in the nonlinear
+ // computes the matrix $N(\cdot,\cdot)$, which appears in the nonlinear
// term in the Jacobian of $F(\cdot)$. Just as <code>compute_nl_term</code>,
// we must allow this function to receive as input an "old" and a "new"
// solution, which we again call $w_{\mathrm{old}}$ and $w_{\mathrm{new}}$
{
make_grid_and_dofs ();
- // To aknowledge the initial condition, we must use the function $u_0(x)$
+ // To acknowledge the initial condition, we must use the function $u_0(x)$
// to compute $U^0$. To this end, below we will create an object of type
// <code>InitialValues</code>; note that when we create this object (which
// is derived from the <code>Function</code> class), we set its internal
// SolutionTransfer class and we have to prepare the solution vectors that
// should be transferred to the new grid (we will lose the old grid once
// we have done the refinement so the transfer has to happen concurrently
- // with refinement). What we definetely need are the current and the old
+ // with refinement). What we definitely need are the current and the old
// temperature (BDF-2 time stepping requires two old solutions). Since the
// SolutionTransfer objects only support to transfer one object per dof
// handler, we need to collect the two temperature solutions in one data
<img src="http://www.dealii.org/images/steps/developer/step-28.grid-1.9.order2.png" alt="">
-We see that the grid of thermal group is much finner than the one of fast group.
+We see that the grid of thermal group is much finer than the one of fast group.
The solutions on these grids are, (Note: flux are normalized with total fission
source equal to 1)
mesh iteration of polynomial order 3 minus 2e-10. We see that h-adaptive calculations
deliver an algebraic convergence. And the higher polynomial order is, the faster mesh
iteration converges. In our problem, we need smaller number of DoFs to achieve same
-accuracy with higher polynoimal order.
+accuracy with higher polynomial order.
// arrays. It takes the number of energy groups as an argument (an throws an
// error if that value is not equal to two, since at presently only data for
// two energy groups is implemented; however, using this, the function
- // remains flexible and extendible into the future). In the member
+ // remains flexible and extendable into the future). In the member
// initialization part at the beginning, it also resizes the arrays to their
// correct sizes.
//
ConstraintMatrix hanging_node_constraints;
- // @sect5{Private member functionss}
+ // @sect5{Private member functions}
//
// There is one private member function in this class. It recursively
// walks over cells of two meshes to compute the cross-group right hand
// groups. First we call <code>get_finest_common_cells</code> to obtain this
// list of pairs of common cells from both meshes. Both cells in a pair may
// not be active but at least one of them is. We then hand each of these
- // cell pairs off to a function tha computes the right hand side terms
+ // cell pairs off to a function that computes the right hand side terms
// recursively.
//
// Note that ingroup coupling is handled already before, so we exit the
will pass through the boundary as if it wasn't here whereas the remaining
fraction of the wave will be reflected back into the domain.
-If we are willing to accept this as a suffient approximation to an absorbing boundary we finally arrive at the following problem for $u$:
+If we are willing to accept this as a sufficient approximation to an absorbing boundary we finally arrive at the following problem for $u$:
@f{eqnarray*}
-\omega^2 u - c^2\Delta u &=& 0, \qquad x\in\Omega,\\
c (n\cdot\nabla u) + i\,\omega\,u &=&0, \qquad x\in\Gamma_2,\\
is ${\cal O}(BN)$. This also explains why the program does run in 3d
as well (after changing the dimension on the
<code>UltrasoundProblem</code> object), but scales very badly and
-takes extraordinate patience before it finishes solving the linear
+takes extraordinary patience before it finishes solving the linear
system on a mesh with appreciable resolution, even though all the
other parts of the program scale very nicely.
{}
- // The actual prostprocessing happens in the following function. Its inputs
+ // The actual postprocessing happens in the following function. Its inputs
// are a vector representing values of the function (which is here
// vector-valued) representing the data vector given to
// DataOut::add_data_vector, evaluated at all evaluation points where we
// At this point, it is important to keep in mind that we are
// dealing with a finite element system with two
// components. Due to the way we constructed this FESystem,
- // namely as the cartesian product of two scalar finite
+ // namely as the Cartesian product of two scalar finite
// element fields, each shape function has only a single
// nonzero component (they are, in deal.II lingo, @ref
// GlossPrimitive "primitive"). Hence, each shape function
common structure that almost all following tutorial programs will use:
<ul>
<li> <code>make_grid()</code>: This is what one could call a
- <i>pre-processing function</i>. As its name suggests, it sets up the
+ <i>preprocessing function</i>. As its name suggests, it sets up the
object that stores the triangulation. In later examples, it could also
deal with boundary conditions, geometries, etc.
<li> <code>setup_system()</code>: This then is the function in which all the
other data structures are set up that are needed to solve the
problem. In particular, it will initialize the DoFHandler object and
correctly size the various objects that have to do with the linear
- algebra. This function is often separated from the pre-processing
+ algebra. This function is often separated from the preprocessing
function above because, in a time dependent program, it may be called
at least every few time steps whenever the mesh
is adaptively refined (something we will see how to do in step-6). On
- the other hand, setting up the mesh itself in the pre-processing
+ the other hand, setting up the mesh itself in the preprocessing
function above is done only once at the beginning of the program and
is, therefore, separated into its own function.
<li> <code>assemble_system()</code>: This, then is where the contents of the
values for those faces that have zero boundary indicator, and leave
those faces alone that have a different boundary indicator. What
this then does is to impose Dirichlet boundary conditions on the
- former, and homogenous Neumann conditions on the latter (i.e. zero
+ former, and homogeneous Neumann conditions on the latter (i.e. zero
normal derivative of the solution, unless one adds additional terms
to the right hand side of the variational equality that deal with
potentially non-zero Neumann conditions). You will see this if you
usual code found in application programs will need modifications to exploit
the features of meshes which were created using anisotropic
refinement. For those interested in how deal.II evolved, it may be of
-interest that the losening of such invariants required some
+interest that the loosening of such invariants required some
incompatible changes. For example, the library used to have a member
GeometryInfo<dim>::children_per_cell that specified how many children
a cell has once it is refined. For isotropic refinement, this number
consequently been removed. It has now been replaced by
GeometryInfo<dim>::max_children_per_cell which specifies the
<i>maximum</i> number of children a cell can have. How many children a
-refined cell has was previously available as static informtion, but
+refined cell has was previously available as static information, but
now it depends on the actual refinement state of a cell and can be
retrieved using the function call <code>cell-@>n_children()</code>,
a call that works equally well for both isotropic and anisotropic
coarse mesh can decrease performance significantly for real world
problems. However, that is not always the case. Considering boundary layers in
compressible viscous flows, for example, the mesh is always aligned with the
-anisotropic features, thus anisotropic refinement will almost alway increase the
+anisotropic features, thus anisotropic refinement will almost always increase the
efficiency of computations on adapted grids for these cases.
// in the left part of the domain at a velocity that matches the one coming
// in from the right. In the circular part the magnitude of the flow
// velocity is proportional to the distance from the origin. This is a
- // difference to step-12, where the magnitude was 1 evereywhere. the new
+ // difference to step-12, where the magnitude was 1 everywhere. the new
// definition leads to a linear variation of $\beta$ along each given face
// of a cell. On the other hand, the solution $u(x,y)$ is exactly the same
// as before.
anisotropic(anisotropic),
// As beta is a linear function, we can choose the degree of the
// quadrature for which the resulting integration is correct. Thus, we
- // choose to use <code>degree+1</code> gauss points, which enables us to
+ // choose to use <code>degree+1</code> Gauss points, which enables us to
// integrate exactly polynomials of degree <code>2*degree+1</code>, enough
// for all the integrals we will perform in this program.
quadrature (degree+1),
// Case b), we decide that there are finer cells as neighbors
// by asking the face, whether it has children. if so, then
// there must also be finer cells which are children or
- // farther offsprings of our neighbor.
+ // farther offspring of our neighbor.
if (face->has_children())
{
// We need to know, which of the neighbors faces points in
(neighbor->level() < cell->level() &&
neighbor->index() == cell->index())))
{
- // Here we know, that the neigbor is not coarser so we
+ // Here we know, that the neighbor is not coarser so we
// can use the usual @p neighbor_of_neighbor
// function. However, we could also use the more
// general @p neighbor_face_no function.
triangulation.execute_coarsening_and_refinement ();
}
- // Once an error indicator has been evaluated and the cells with largerst
+ // Once an error indicator has been evaluated and the cells with largest
// error are flagged for refinement we want to loop over the flagged cells
- // again to decide whether they need isotropic refinemnt or whether
+ // again to decide whether they need isotropic refinement or whether
// anisotropic refinement is more appropriate. This is the anisotropic jump
// indicator explained in the introduction.
template <int dim>
endc=dof_handler.end();
for (; cell!=endc; ++cell)
- // We only need to consider cells which are flaged for refinement.
+ // We only need to consider cells which are flagged for refinement.
if (cell->refine_flag_set())
{
Point<dim> jump;
std::vector<double> u (fe_v_face.n_quadrature_points);
std::vector<double> u_neighbor (fe_v_face.n_quadrature_points);
- // The four cases of different neighbor relations senn in
+ // The four cases of different neighbor relations seen in
// the assembly routines are repeated much in the same way
// here.
if (face->has_children())
// present subface...
typename DoFHandler<dim>::cell_iterator neighbor_child = cell->neighbor_child_on_subface(face_no,subface_no);
Assert (!neighbor_child->has_children(), ExcInternalError());
- // ... and reinit the respective FEFaceValues und
+ // ... and reinit the respective FEFaceValues and
// FESubFaceValues objects.
fe_v_subface.reinit (cell, face_no, subface_no);
fe_v_face_neighbor.reinit (neighbor_child, neighbor2);
fe_v_subface.get_function_values(solution2, u);
fe_v_face_neighbor.get_function_values(solution2, u_neighbor);
// as well as the quadrature weights, multiplied by
- // the jacobian determinant.
+ // the Jacobian determinant.
const std::vector<double> &JxW = fe_v_subface.get_JxW_values ();
// Now we loop over all quadrature points
for (unsigned int x=0; x<fe_v_subface.n_quadrature_points; ++x)
// first coordinate direction on the unit cell,
// the second two faces are orthogonal to the
// second coordinate direction and so on, so we
- // accumulate these values ito vectors with
+ // accumulate these values into vectors with
// <code>dim</code> components.
jump[face_no/2]+=std::fabs(u[x]-u_neighbor[x])*JxW[x];
// We also sum up the scaled weights to obtain
// Now we loop over the <code>dim</code> coordinate directions of
// the unit cell and compare the average jump over the faces
- // orthogional to that direction with the average jumnps over faces
- // orthogonal to the remining direction(s). If the first is larger
+ // orthogonal to that direction with the average jumps over faces
+ // orthogonal to the remaining direction(s). If the first is larger
// than the latter by a given factor, we refine only along hat
// axis. Otherwise we leave the refinement flag unchanged, resulting
// in isotropic refinement.
melting point, $\eta$ may change by orders of magnitude over the typical range
of temperatures.
-We note that the Stokes equation above could be non-dimensionalized by
+We note that the Stokes equation above could be nondimensionalized by
introducing the <a target="_top"
href="http://en.wikipedia.org/wiki/Rayleigh_number">Rayleigh
number</a> $\mathrm{Ra}=\frac{\|g\| \beta \rho}{\eta \kappa} \delta T L^3$ using a
numerical viscosity $\nu(T)$ is temperature-dependent, so the equation
we want to solve is nonlinear in <i>T</i> — not what one desires from a
simple method to stabilize an equation, and even less so if we realize
-that $\nu(T)$ is non-differentiable in <i>T</i>. However, there is no
+that $\nu(T)$ is nondifferentiable in <i>T</i>. However, there is no
reason to despair: we still have to discretize in time and we can
treat the term explicitly.
temperature at the current time in the same way as we have discussed
the advection works in the section on time stepping.
-The form for non-uniform time steps that we will have to use in
+The form for nonuniform time steps that we will have to use in
reality is a bit more complicated (which is why we showed the simpler
form above first) and reads:
@f{eqnarray*}
introduction of step-22 and fixes one component of the velocity
while allowing flow to be %parallel to the boundary. There remain
<code>dim-1</code> components to be fixed, namely the tangential components of
-the normal stress; for these, we choose homogenous conditions which means that
+the normal stress; for these, we choose homogeneous conditions which means that
we do not have to anything special. Initial conditions are only necessary for
the temperature field, and we choose it to be constant zero.
close to the bottom of the domain, imitating three candles. The fluid located
at these sources, initially at rest, is then heated up and as the temperature
rises gains buoyancy, rising up; more fluid is dragged up and through the
-sources, leading to three hote plumes that rise up until they are captured by
+sources, leading to three hot plumes that rise up until they are captured by
the recirculation of fluid that sinks down on the outside, replacing the air
that rises due to heating.
small enough so that the distance transport advects in each time step
is no longer than one <i>grid point</i> away (which for $Q_1$ elements
is $h_K$, but for $Q_2$ elements is $h_K/2$). It turns out that $\beta$
-needs to be sligthly larger for obtaining stable results also late in
+needs to be slightly larger for obtaining stable results also late in
the simulation at times larger than 60, so we actually choose it as
$\beta = 0.034$ in the code.
// same as in step-22. Regarding the details, though, there are some
// differences.
- // The first thing is that we don't set any non-homogenous boundary
+ // The first thing is that we don't set any nonhomogeneous boundary
// conditions on the velocity, since as is explained in the introduction we
// will use no-flux conditions $\mathbf{n}\cdot\mathbf{u}=0$. So what is
// left are <code>dim-1</code> conditions for the tangential part of the
// normal component of the stress tensor, $\textbf{n} \cdot [p \textbf{1} -
- // \eta\varepsilon(\textbf{u})]$; we assume homogenous values for these
+ // \eta\varepsilon(\textbf{u})]$; we assume homogeneous values for these
// components, i.e. a natural boundary condition that requires no specific
// action (it appears as a zero term in the right hand side of the weak
// form).
//
// Secondly, we catch any exceptions that the solver may have thrown. The
// reason is as follows: When debugging a program like this one
- // occasionally makes a mistake of passing an indefinite or non-symmetric
+ // occasionally makes a mistake of passing an indefinite or nonsymmetric
// matrix or preconditioner to the current class. The solver will, in that
// case, not converge and throw a run-time exception. If not caught here
// it will propagate up the call stack and may end up in
// detail in the results section of this tutorial program. The second is the
// exponent $\alpha$; $\alpha=1$ appears to work fine for the current
// program, even though some additional benefit might be expected from
- // chosing $\alpha = 2$. Finally, there is one thing that requires special
+ // choosing $\alpha = 2$. Finally, there is one thing that requires special
// casing: In the first time step, the velocity equals zero, and the formula
// for $\nu|_K$ is not defined. In that case, we return $\nu|_K=5\cdot 10^3
// \cdot h_K$, a choice admittedly more motivated by heuristics than
// dofs. Consequently, there will be no data written at positions of
// constrained degrees of freedom, so we can let the
// DoFTools::make_sparsity_pattern function omit these entries by setting
- // the last boolean flag to <code>false</code>. Once the sparsity pattern
+ // the last Boolean flag to <code>false</code>. Once the sparsity pattern
// is ready, we can use it to initialize the Trilinos matrices. Since the
// Trilinos matrices store the sparsity pattern internally, there is no
// need to keep the sparsity pattern around after the initialization of
// Next, we set some more options of the AMG preconditioner. In
// particular, we need to tell the AMG setup that we use quadratic basis
// functions for the velocity matrix (this implies more nonzero elements
- // in the matrix, so that a more rubust algorithm needs to be chosen
+ // in the matrix, so that a more robust algorithm needs to be chosen
// internally). Moreover, we want to be able to control how the coarsening
// structure is build up. The way the Trilinos smoothed aggregation AMG
// does this is to look which matrix entries are of similar size as the
// diagonal entry in order to algebraically build a coarse-grid
// structure. By setting the parameter <code>aggregation_threshold</code>
- // to 0.02, we specify that all entries that are more than two precent of
+ // to 0.02, we specify that all entries that are more than two percent of
// size of some diagonal pivots in that row should form one coarse grid
- // point. This parameter is rather ad-hoc, and some fine-tuning of it can
+ // point. This parameter is rather ad hoc, and some fine-tuning of it can
// influence the performance of the preconditioner. As a rule of thumb,
// larger values of <code>aggregation_threshold</code> will decrease the
// number of iterations, but increase the costs per iteration. A look at
//
// Finally, we also initialize the preconditioner for the inversion of the
// pressure mass matrix. This matrix is symmetric and well-behaved, so we
- // can chose a simple preconditioner. We stick with an incomple Cholesky
+ // can chose a simple preconditioner. We stick with an incomplete Cholesky
// (IC) factorization preconditioner, which is designed for symmetric
// matrices. We could have also chosen an SSOR preconditioner with
// relaxation factor around 1.2, but IC is cheaper for our example. We
// Next comes the declaration of vectors to hold the old and older
// solution values (as a notation for time levels <i>n-1</i> and
// <i>n-2</i>, respectively) and gradients at quadrature points of the
- // current cell. We also declarate an object to hold the temperature right
- // hande side values (<code>gamma_values</code>), and we again use
+ // current cell. We also declare an object to hold the temperature right
+ // hand side values (<code>gamma_values</code>), and we again use
// shortcuts for the temperature basis functions. Eventually, we need to
// find the temperature extrema and the diameter of the computational
// domain which will be used for the definition of the stabilization
// according to the discussion in the introduction using the dedicated
// function. With that at hand, we can get into the loop over
// quadrature points and local rhs vector components. The terms here
- // are quite lenghty, but their definition follows the time-discrete
+ // are quite lengthy, but their definition follows the time-discrete
// system developed in the introduction of this program. The BDF-2
// scheme needs one more term from the old time step (and involves
// more complicated factors) than the backward Euler scheme that is
// There is a snatch here. The formula contains a division by the maximum
// value of the velocity. However, at the start of the computation, we
// have a constant temperature field (we start with a constant
- // temperature, and it will be non-constant only after the first time step
+ // temperature, and it will be nonconstant only after the first time step
// during which the source acts). Constant temperature means that no
// buoyancy acts, and so the velocity is zero. Dividing by it will not
// likely lead to anything good.
// SolutionTransfer class and we have to prepare the solution vectors that
// should be transferred to the new grid (we will lose the old grid once
// we have done the refinement so the transfer has to happen concurrently
- // with refinement). What we definetely need are the current and the old
+ // with refinement). What we definitely need are the current and the old
// temperature (BDF-2 time stepping requires two old solutions). Since the
// SolutionTransfer objects only support to transfer one object per dof
// handler, we need to collect the two temperature solutions in one data
reason back there is that everything was nicely balanced: velocities
were on the order of one, the pressure likewise, the viscosity was
one, and the domain had a diameter of $\sqrt{2}$. As a result, while
-non-sensical, nothing bad happened. On the other hand, as we will explain
+nonsensical, nothing bad happened. On the other hand, as we will explain
below, things here will not be that simply scaled: $\eta$ will be around
$10^{21}$, velocities on the order of $10^{-8}$, pressure around $10^8$, and
the diameter of the domain is $10^7$. In other words, the order of magnitude
matrix and preconditioner. Because it is annoying and error prone, we will
recover the unscaled pressure immediately following the solution of the linear
system, i.e., the solution vector's pressure component will immediately be
-un-scaled to retrieve the physical pressure. Since the solver uses the fact that
+unscaled to retrieve the physical pressure. Since the solver uses the fact that
we can use a good initial guess by extrapolating the previous solutions, we
also have to scale the pressure immediately <i>before</i> solving.
with the unreasonably coarse mesh chosen for these experiments, around
time $t=10^{15}$ seconds the plumes of hot material that have been
rising towards the cold outer boundary and have then spread sideways
-are starting to get close to each other, squezzing out the cold
-material inbetween. This creates a layer of cells into which fluids
+are starting to get close to each other, squeezing out the cold
+material in-between. This creates a layer of cells into which fluids
flows from two opposite sides and flows out toward a third, apparently
a scenario that then produce these instabilities without sufficient
stabilization. Second: In step-31, we used
physics would have to be implemented, and several other aspects are
currently missing from this program as well. We will come back to this
issue in the results section again, but state for now that providing a
-realistic description is a goal of the <i>Apect</i> code in
+realistic description is a goal of the <i>Aspect</i> code in
development at the time of writing this.
As a reminder, let us again state the equations we want to solve are these:
entry</a>). In Kelvin, we therefore choose $T_0=(4000+273)\text{K}$,
$T_1=(500+273)\text{K}$ as boundary conditions at the inner and outer edge.
- In addition to this, we also have to specifiy some initial conditions for
+ In addition to this, we also have to specify some initial conditions for
the temperature field. The real temperature field of the earth is quite
complicated as a consequence of the convection that has been going on for
more than four billion years -- in fact, it is the properties of this
@f]
- There are two problems with this, however: (i) The Earth is not homogenous,
+ There are two problems with this, however: (i) The Earth is not homogeneous,
i.e. the density $\rho$ depends on $\mathbf x$; in fact it is not even a
function that only depends on the radius $r=\|\mathbf x\|$. In reality, gravity therefore
does not always decrease as we get deeper: because the earth core is so much
changes as temperature and pressure change: depending on temperature and
pressure, different crystal configurations are thermodynamically favored
over others, even if the chemical composition of the mantle were
- homogenous. For example, the common mantle material MgSiO<sub>3</sub> exists
+ homogeneous. For example, the common mantle material MgSiO<sub>3</sub> exists
in its <a target="_top"
href="http://en.wikipedia.org/wiki/Perovskite_(structure)">perovskite
structure</a> throughout most of the mantle, but in the lower mantle the
program just because one wants to play with a single parameter (think, for
example, of parameter studies determining the best values of the
stabilization constants discussed above), in particular given that it takes
- a non-trivial amount of time to re-compile programs of the current size. To
+ a nontrivial amount of time to re-compile programs of the current size. To
just give an overview of the kinds of parameters we have moved from fixed
values into the input file, here is a listing of a typical
<code>\step-32.prm</code> file:
trivial. An additional complexity is introduced by the fact that one may
want to change the number of processors between runs, for example because
one may wish to continue computing on a mesh that is finer than the one used
- to pre-compute a starting temperature field at an intermediate time.
+ to precompute a starting temperature field at an intermediate time.
<li> <b>Predictive postprocessing:</b> The point of computations like this is
not simply to solve the equations. Rather, it is typically the exploration
// In the following namespace, we define the various pieces of equation data
// that describe the problem. This corresponds to the various aspects of
- // making the problem at least slightly realistc and that were exhaustively
+ // making the problem at least slightly realistic and that were exhaustively
// discussed in the description of the testcase in the introduction.
//
// We start with a few coefficients that have constant values (the comment
// block. The three code blocks of the <code>vmult</code> function implement
// the multiplications with the three blocks of this preconditioner matrix
// and should be self explanatory if you have read through step-31 or the
- // discussion of compositing solvers in step-20.
+ // discussion of composing solvers in step-20.
namespace LinearSolvers
{
template <class PreconditionerA, class PreconditionerMp>
}
// Following this, we can compute constraints for the solution vectors,
- // including hanging node constraints and homogenous and inhomogenous
+ // including hanging node constraints and homogeneous and inhomogeneous
// boundary values for the Stokes and temperature fields. Note that as for
// everything else, the constraint objects can not hold <i>all</i>
// constraints on every processor. Rather, each processor needs to store
// dof on a cell, the joint finite element knows to which equation component
// (velocity component, pressure, or temperature) it belongs – that's the
// information we need! So we step through all cells (with iterators into
- // all three DoFHandlers moving in synch), and for each joint cell dof, we
+ // all three DoFHandlers moving in sync), and for each joint cell dof, we
// read out that component using the FiniteElement::system_to_base_index
// function (see there for a description of what the various parts of its
// return value contain). We also need to keep track whether we're on a
The output are the derivatives $\frac{\partial c(a,b)}{\partial a},
\frac{\partial c(a,b)}{\partial b}$ of $c(a,b)=2a+\cos(ab)$ at $a=1,b=2$.
-It should be noted that Sacado provides more auto-differentation capabilities than the small subset
+It should be noted that Sacado provides more auto-differentiation capabilities than the small subset
used in this program. However, understanding the example above is
enough to understand the use of Sacado in this Euler flow program.
within the example.
<h3> Adaptivity </h3>
-The example uses an ad-hoc refinement indicator that shows some usefulness in shock-type problems, and
+The example uses an ad hoc refinement indicator that shows some usefulness in shock-type problems, and
in the downhill flow example included. We refine according to the squared gradient of the density.
Hanging nodes are handled by computing the numerical flux across cells that are of differing
refinement levels, rather than using the ConstraintMatrix class as in
{
// We prescribe the velocity (we are dealing with a particular
// component here so that the average of the velocities is
- // orthogonal to the surface normal. This creates sensitivies of
+ // orthogonal to the surface normal. This creates sensitivities of
// across the velocity components.
Sacado::Fad::DFad<double> vdotn = 0;
for (unsigned int d = 0; d < dim; d++)
// This is the only function worth commenting on. When generating graphical
// output, the DataOut and related classes will call this function on each
- // cell, with values, gradients, hessians, and normal vectors (in case we're
+ // cell, with values, gradients, Hessians, and normal vectors (in case we're
// working on faces) at each quadrature point. Note that the data at each
// quadrature point is itself vector-valued, namely the conserved
// variables. What we're going to do here is to compute the quantities we're
// interested in at each quadrature point. Note that for this we can ignore
- // the hessians ("dduh") and normal vectors; to avoid compiler warnings
+ // the Hessians ("dduh") and normal vectors; to avoid compiler warnings
// about unused variables, we comment out their names.
template <int dim>
void
// we allow up to <code>max_n_boundaries</code> boundary indicators to be
// used in the input file, and each of these boundary indicators can be
// associated with an inflow, outflow, or pressure boundary condition,
- // with inhomogenous boundary conditions being specified for each
+ // with homogeneous boundary conditions being specified for each
// component and each boundary indicator separately.
//
// The data structure used to store the boundary indicators is a bit
// This final set of member variables (except for the object holding all
// run-time parameters at the very bottom and a screen output stream that
// only prints something if verbose output has been requested) deals with
- // the inteface we have in this program to the Trilinos library that
+ // the interface we have in this program to the Trilinos library that
// provides us with linear solvers. Similarly to including PETSc matrices
// in step-17, step-18, and step-19, all we need to do is to create a
// Trilinos sparse matrix instead of the standard deal.II class. The
// whereas all the other ones remain dependent functions. These are
// precisely the local degrees of freedom just extracted. All calculations
// that reference them (either directly or indirectly) will accumulate
- // sensitivies with respect to these variables.
+ // sensitivities with respect to these variables.
//
// In order to mark the variables as independent, the following does the
// trick, marking <code>independent_local_dof_values[i]</code> as the
// (\mathbf{z}_i)_{\text{component\_i}})_K$, where integrals are
// understood to be evaluated through summation over quadrature points.
//
- // We initialy sum all contributions of the residual in the positive
+ // We initially sum all contributions of the residual in the positive
// sense, so that we don't need to negative the Jacobian entries. Then,
// when we sum into the <code>right_hand_side</code> vector, we negate
// this residual.
// Now assemble the face term in exactly the same way as for the cell
// contributions in the previous function. The only difference is that if
// this is an internal face, we also have to take into account the
- // sensitivies of the residual contributions to the degrees of freedom on
+ // sensitivities of the residual contributions to the degrees of freedom on
// the neighboring cell:
std::vector<double> residual_derivatives (dofs_per_cell);
for (unsigned int i=0; i<fe_v.dofs_per_cell; ++i)
// Note that the number of the output file is determined by keeping a
// counter in the form of a static variable that is set to zero the first
// time we come to this function and is incremented by one at the end of
- // each invokation.
+ // each invocation.
template <int dim>
void ConservationLaw<dim>::output_results () const
{
\mathbf{\tilde{v}}=\nabla\phi,
\f]
and so the second part of Euler equations above can be rewritten
-as the homogenous Laplace equation for the unknown $\phi$:
+as the homogeneous Laplace equation for the unknown $\phi$:
\f{align*}
\label{laplace}
\Delta\phi &= 0 \qquad &&\text{in}\ \mathbb{R}^n\backslash\Omega,
If this were a production code rather than a demonstration of principles,
there are techniques that are available to not store full matrices but instead
store only those entries that are large and/or relevant. In the literature on
-boundary element methods, a plethory of methods is available that allows to
+boundary element methods, a plethora of methods is available that allows to
determine which elements are important and which are not, leading to a
significantly sparser representation of these matrices that also facilitates
rapid evaluations of the scalar product between vectors and matrices. This not
<img src="http://www.dealii.org/images/steps/developer/step-34_2d.png" alt="">
while in three dimensions we show first the potential on the surface,
-together with a contur plot,
+together with a contour plot,
<img src="http://www.dealii.org/images/steps/developer/step-34_3d.png" alt="">
// @sect4{BEMProblem::BEMProblem and BEMProblem::read_parameters}
- // The constructor initializes the variuous object in much the same way as
+ // The constructor initializes the various object in much the same way as
// done in the finite element programs such as step-4 or step-6. The only
// new ingredient here is the ParsedFunction object, which needs, at
// construction time, the specification of the number of components.
// the current cell, a deviation from the usual finite element
// integrals). The integral that we need to perform is singular if one
// of the local degrees of freedom is the same as the support point
- // $i$. A the beginning of the loop we therefore check wether this is
+ // $i$. A the beginning of the loop we therefore check whether this is
// the case, and we store which one is the singular index:
for (unsigned int i=0; i<dh.n_dofs() ; ++i)
{
// Singular integration requires a careful selection of the quadrature
// rules. In particular the deal.II library provides quadrature rules which
- // are taylored for logarithmic singularities (QGaussLog, QGaussLogR), as
+ // are tailored for logarithmic singularities (QGaussLog, QGaussLogR), as
// well as for 1/R singularities (QGaussOneOverR).
//
// Singular integration is typically obtained by constructing weighted
// consideration also the second part of the integral.
//
// A similar reasoning should be done in the three dimensional case, since
- // the singular quadrature is taylored on the inverse of the radius $r$ in
+ // the singular quadrature is tailored on the inverse of the radius $r$ in
// the reference cell, while our singular function lives in real space,
// however in the three dimensional case everything is simpler because the
// singularity scales linearly with the determinant of the
// The implementation of a diffusion step. Note that the expensive operation
// is the diffusion solve at the end of the function, which we have to do
- // once for each velocity component. To accellerate things a bit, we allow
+ // once for each velocity component. To accelerate things a bit, we allow
// to do this in %parallel, using the Threads::new_task function which makes
// sure that the <code>dim</code> solves are all taken care of and are
// scheduled to available processors: if your machine has more than one
Schrödinger equation which models the motion of particles in an
external potential $V(\mathbf x)$. The particle is described by a wave
function $\Psi(\mathbf x)$ that satisfies a relation of the
-(non-dimensionalized) form
+(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.
the various eigenfunctions we want to compute. Moreover, the
<code>main</code> function, which has the top-level control over
everything here, initializes and finalizes the interface to SLEPc and
-PETSc simultaneouly via <code>SlepcInitialize</code> and
+PETSc simultaneously via <code>SlepcInitialize</code> and
<code>SlepFinalize</code>.</li>
<li>We use PETSc matrices and vectors as in step-17 and
machine architecture. However, for larger grids and with a larger
number of degrees-of-freedom, we may want to run our application on
parallel architectures. A parallel implementation of the above code
-can be particularily useful here since the generalized eigenspectrum
+can be particularly useful here since the generalized eigenspectrum
problem is somewhat more expensive to solve than the standard problems
considered in most of the earlier tutorials. Fortunately, modifying the above
program to be MPI compliant is a relatively straightforward
This example shows how to implement a matrix-free method, that is, a
method that does not explicitly store the matrix elements, for a
second-order Poisson equation with variable coefficients on a
-hypercube. The eliptic equation will be solved with a multigrid
+hypercube. The elliptic equation will be solved with a multigrid
method.
The major motivation for matrix-free methods is the fact that today
for (unsigned int q=0; q<n_q_points; ++q)
{
- // apply the tranpose inverse Jacobian of the mapping
+ // apply the transposed inverse Jacobian of the mapping
Tensor<1,dim> temp;
for (unsigned int d=0; d<dim; ++d)
temp[d] = temp_vector(q*dim+d);
linear elements, but faster for all higher order elements thanks to the
reduced complexity due to the tensor product structure and due to less memory
transfer during computations. The impact of reduced memory transfer is
-particularly beneficial when working on a multi-core processor where several
+particularly beneficial when working on a multicore processor where several
processing units share access to memory. In that case, an algorithm which is
computation bound will show almost perfect parallel speedup, whereas an
algorithm that is bound by memory transfer might not achieve similar speedup
// would have to have two classes like the current one for each of the
// operators (maybe with a common base class). However, in that case, we
// would not store a MatrixFree object in this class to avoid doing the
- // expensive work of pre-computing everything MatrixFree stores
+ // expensive work of precomputing everything MatrixFree stores
// twice. Rather, we would keep this object in the main class and simply
// store a reference.
//
// @sect4{Initialization}
- // Once we have created the multi-grid dof_handler and the constraints, we
- // can call the reinit function for each level of the multi-grid routine
+ // Once we have created the multigrid dof_handler and the constraints, we
+ // can call the reinit function for each level of the multigrid routine
// (and the active cells). The main purpose of the reinit function is to
// setup the <code> MatrixFree </code> instance for the problem. Also, the
// coefficient is evaluated. For this, we need to activate the update flag
MappingQ that will automatically obtain the boundary description from the
Triangulation. The mapping object will then be passed to the appropriate
functions, and we will get a boundary description for half circles or half
-spheres that is pre-defined in the library.
+spheres that is predefined in the library.
The rest of the program follows closely step-4 and, as far as computing the
error, step-7. Some aspects of this program, in particular the use of two
Note that the only essential addition has been the three lines marked with
asterisks. It is worth pointing out one other thing here, though: because we
-un-attach the manifold description from the surface mesh, whenever we use a
+detach the manifold description from the surface mesh, whenever we use a
mapping object in the rest of the program, it has no curves boundary
description to go on any more. Rather, it will have to use the implicit,
StraightBoundary class that is used on all parts of the boundary not
-explicitly assigned a different mannifold object. Consequently, whether we use
+explicitly assigned a different manifold object. Consequently, whether we use
MappingQ(2), MappingQ(15) or MappingQ1, each cell of our mesh will be mapped
using a bilinear approximation.
// boundary of the domain -- this is what we did in step-11, for
// example. However, here we have a curved domain, not just a curved
// boundary, and while we can approximate it with bilinearly mapped cells,
- // it is really only prodent to use a higher order mapping for all
+ // it is really only prudent to use a higher order mapping for all
// cells. Consequently, this class has a member variable of type MappingQ;
// we will choose the polynomial degree of the mapping equal to the
// polynomial degree of the finite element used in the computations to
<h1>Results</h1>
<h2>Logfile output</h2>
-First, the program produces the usual logfile here stored in <tt>deallog</tt>. It reads (with ommission of intermediate steps)
+First, the program produces the usual logfile here stored in <tt>deallog</tt>. It reads (with omission of intermediate steps)
@code
DEAL::Element: FE_DGQ<2>(3)
// ourselves. Similarly, we implement Nitsche boundary conditions and the
// interior penalty fluxes between cells.
//
- // The boundary und flux terms need a penalty parameter, which should be
+ // The boundary and flux terms need a penalty parameter, which should be
// adjusted to the cell size and the polynomial degree. A safe choice of
// this parameter for constant coefficients can be found in
// LocalIntegrators::Laplace::compute_penalty() and we use this below.
};
// Here we have the integration on cells. There is currently no good
- // interfce in MeshWorker that would allow us to access values of regular
+ // interface in MeshWorker that would allow us to access values of regular
// functions in the quadrature points. Thus, we have to create the vectors
// for the exact function's values and gradients inside the cell
// integrator. After that, everything is as before and we just add up the
mg_matrix_dg_down.clear();
// It is important to update the sparsity patterns after <tt>clear()</tt>
// was called for the level matrices, since the matrices lock the sparsity
- // pattern through the Smartpointer ans Subscriptor mechanism.
+ // pattern through the Smartpointer and Subscriptor mechanism.
mg_sparsity.resize(0, n_levels-1);
mg_sparsity_dg_interface.resize(0, n_levels-1);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
- // Since this assembler alows us to fill several vectors, the interface is
+ // Since this assembler allows us to fill several vectors, the interface is
// a little more complicated as above. The pointers to the vectors have to
// be stored in a NamedData object. While this seems to cause two extra
// lines of code here, it actually comes handy in more complex
const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1;
info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points);
- // but now we need to notify the info box of the finite element functio we
+ // but now we need to notify the info box of the finite element function we
// want to evaluate in the quadrature points. First, we create a NamedData
// object with this vector, which is the solution we just computed.
NamedData<Vector<double>* > solution_data;
// Then, we tell the Meshworker::VectorSelector for cells, that we need
// the second derivatives of this solution (to compute the
- // Laplacian). Therefore, the boolean arguments selecting function values
+ // Laplacian). Therefore, the Boolean arguments selecting function values
// and first derivatives a false, only the last one selecting second
// derivatives is true.
info_box.cell_selector.add("solution", false, false, true);
and 3d!). We would have to write the same function twice, once
working on <code>Triangulation2d</code> and once working with a
<code>Triangulation3d</code>. This is an unnecessary obstacle in
-programming and leads to a nuisance to keep the two function in synch
+programming and leads to a nuisance to keep the two function in sync
(at best) or difficult to find errors if the two versions get out of
sync (at worst; this would probably the more common case).
what is happening on the inside. The picture on the right is an attempt at
visualizing the interior as well, by showing surfaces where the solution has
constant values (as indicated by the legend at the top left). Isosurface
-pictures look best if one makes the individual surfaces slightly transparaent
+pictures look best if one makes the individual surfaces slightly transparent
so that it is possible to see through them and see what's behind.
<TABLE WIDTH="60%" ALIGN="center">
// formulas above right away.
//
// The last thing to note is that a <code>Point@<dim@></code> denotes a point
-// in dim-dimensionsal space, and its individual components (i.e. $x$, $y$,
+// in dim-dimensional space, and its individual components (i.e. $x$, $y$,
// ... coordinates) can be accessed using the () operator (in fact, the []
// operator will work just as well) with indices starting at zero as usual in
// C and C++.
// <code>data_out.write_gnuplot</code> call by
// <code>data_out.write_vtk</code>.
//
-// Since the program will run both 2d and 3d versions of the laplace solver,
+// Since the program will run both 2d and 3d versions of the Laplace solver,
// we use the dimension in the filename to generate distinct filenames for
// each run (in a better program, one would check whether <code>dim</code> can
// have other values than 2 or 3, but we neglect this here for the sake of
this problem, the
limiting factor is that the program uses the BoomerAMG algebraic
multigrid method from the <a
-href="http://acts.nersc.gov/hypre/" target="_top">hypre package</a> as
+href="http://acts.nersc.gov/hypre/" target="_top">Hypre package</a> as
a preconditioner, which unfortunately uses signed 32-bit integers to
index the elements of a %distributed matrix. This limits the size of
problems to $2^31-1=2,147,483,647$ degrees of freedom. From the graphs
dof_handler.distribute_dofs (fe);
- // The next two lines extract some informatino we will need later on,
+ // The next two lines extract some information we will need later on,
// namely two index sets that provide information about which degrees of
// freedom are owned by the current processor (this information will be
// used to initialize solution and right hand side vectors, and the system
// that stores, for each cell, the subdomain the cell belongs to. This is
// slightly tricky, because of course not every processor knows about every
// cell. The vector we attach therefore has an entry for every cell that the
- // current processor has in its mesh (locally owned onces, ghost cells, and
+ // current processor has in its mesh (locally owned ones, ghost cells, and
// artificial cells), but the DataOut class will ignore all entries that
// correspond to cells that are not owned by the current processor. As a
// consequence, it doesn't actually matter what values we write into these
// The final function, <code>main()</code>, again has the same structure as in
// all other programs, in particular step-6. Like in the other programs that
-// use PETSc, we have to inialize and finalize PETSc, which is done using the
+// use PETSc, we have to initialize and finalize PETSc, which is done using the
// helper object MPI_InitFinalize.
//
// Note how we enclose the use the use of the LaplaceProblem class in a pair
// of braces. This makes sure that all member variables of the object are
-// destroyed by the time we destroy the mpi_intialization object. Not doing
+// destroyed by the time we destroy the mpi_initialization object. Not doing
// this will lead to strange and hard to debug errors when
// <code>PetscFinalize</code> first deletes all PETSc vectors that are still
// around, and the destructor of the LaplaceProblem class then tries to delete
stair step obstacle against which gravity pushes the membrane.
This problem is typically called the "obstacle problem" (see also <a
-href="http://en.wikipedia.org/wiki/Obstacle_problem">this wikipedia article</a>), and it results in a
+href="http://en.wikipedia.org/wiki/Obstacle_problem">this Wikipedia article</a>), and it results in a
variational inequality, rather than a variational equation when put into the
weak form. We will below derive it from the classical formulation, but before we
go on to discuss the mathematics let us show how the solution of the problem we
@f{equation*}
E(u)\leq E(v)\quad \forall v\in G,
@f}
-with the convex set of admissble displacements:
+with the convex set of admissible displacements:
@f{equation*}
G:=\lbrace v\in V: v\geq g \text{ a.e. in } \Omega\rbrace,\quad V:=H^1_0(\Omega).
@f}
On the other hand, if $u=g$ then $G$ only allows test functions $v$ so that in fact
$v-u\ge 0$. This means that we can't test the equation with both $v-u$ and
$-(v-u)$ as above, and so we can no longer conclude that the two sides are in
-fact equal. Thus, this mimicks the way we have discussed the complementarity
+fact equal. Thus, this mimics the way we have discussed the complementarity
condition above.
large enough; for example there is no convergence for $c = 1$ using the
current program if we use 7 global refinements.)
-After some headscratching one can then convince oneself that the inequalities
+After some head-scratching one can then convince oneself that the inequalities
above can equivalently be rewritten as
@f{equation*}
C([BU]_i,\Lambda_i) = 0, \qquad \forall i\in{\cal S}.
// function given above and in the introduction.
//
// If we decide that the DoF should be part of the active set, we
- // add its index to the active set, introduce a nonhomogeneous
+ // add its index to the active set, introduce an inhomogeneous
// equality constraint in the ConstraintMatrix object, and reset the
// solution value to the height of the obstacle. Finally, the
// residual of the non-contact part of the system serves as an
// @sect4{ObstacleProblem::output_results}
// We use the vtk-format for the output. The file contains the displacement
- // and a numerical represenation of the active set. The function looks
+ // and a numerical representation of the active set. The function looks
// standard but note that we can add an IndexSet object to the DataOut
// object in exactly the same way as a regular solution vector: it is simply
// interpreted as a function that is either zero (when a degree of freedom
The continuous and convex function $\mathcal{F}$ denotes the von Mises flow function
$$\mathcal{F}(\tau) = \vert\tau^D\vert - \sigma_0¸\quad\text{with}\quad \tau^D
= \tau - \dfrac{1}{3}tr(\tau)I,$$
-$\sigma_0$ as yield stress and $\vert .\vert$ as the frobenius norm. If there
+$\sigma_0$ as yield stress and $\vert .\vert$ as the Frobenius norm. If there
are no plastic deformations in a particular point - that is $\lambda=0$ - this yields $\vert\sigma^D\vert <
\sigma_0$ and otherwise if $\lambda > 0$ it follows that $\vert\sigma^D\vert = \sigma_0$.
That means if the stress is smaller than the yield stress there are only elastic
deformations in that point.\\
To consider it the other way around if the deviator stress $\sigma^D$ is in a
norm bigger than the yield stress then $\sigma^D$ has to be projected back to the yield surface and there are plastic deformations which means $\lambda$
-would be positiv for that particular point. We refer that the stresses are
-computed by Hooke's law for isotorpic materials. You can find the description at the end of section 3. Else if the norm of the deviator stress tensor is smaller or equal the yield stress then $\lambda$ is zero and there are no plastic deformations in
+would be positive for that particular point. We refer that the stresses are
+computed by Hooke's law for isotropic materials. You can find the description at the end of section 3. Else if the norm of the deviator stress tensor is smaller or equal the yield stress then $\lambda$ is zero and there are no plastic deformations in
that point.\\
There the index $D$ denotes the deviator part of for example the stress where
$tr(.)$ is the trace of a tensor. The definition shows an additive decomposition
density $f$ which we will neglect in our example.
The boundary of $\Omega$ separates as follows $\Gamma=\Gamma_D\bigcup\Gamma_C$ and $\Gamma_D\bigcap\Gamma_C=\emptyset$.
At the boundary $\Gamma_D$ we have zero Dirichlet conditions. $\Gamma_C$ denotes the potential contact boundary.\\
-The last two lines decribe the so-called Signorini contact conditions. If there is no contact the normal stress
+The last two lines describe the so-called Signorini contact conditions. If there is no contact the normal stress
$$ \sigma_n = \sigma n\cdot n$$
is zero with the outward normal $n$. If there is contact ($u_n = g$) the tangential stress $\sigma_t = \sigma\cdot n - \sigma_n n$
vanishes, because we consider a frictionless situation and the normal stress is
L^2(\Omega,\mathbb{R}^{\textrm{dim}\times\textrm{dim}}_{\textrm{sym}}):\textrm{div}(\tau)\in L^2(\Omega,\mathbb{R}^{\textrm{dim}})\rbrace$$ and
$$\Pi \Sigma:=\lbrace \tau\in \Sigma, \mathcal{F}(\tau)\leq 0\rbrace$$
as the set of admissible stresses which is defined
-by a continious, convex flow function $\mathcal{F}$.
+by a continuous, convex flow function $\mathcal{F}$.
With the goal of deriving the dual formulation of the minimisation
-problem, we define a lagrange function:
+problem, we define a Lagrange function:
$$L(\tau,\varphi) := E(\tau) + (\varphi, \textrm{div}(\tau)),\quad \lbrace\tau,\varphi\rbrace\in\Pi W^{\textrm{div}}\times V^+$$
with
$$V^+ := \lbrace u\in V: u_n\leq g \text{ on } \Gamma_C \rbrace$$
$$\tau^D := C\varepsilon^D(u^i).$$
Remark that $a(.;.)$ is not differentiable in the common sense but it is
slantly differentiable like the function for the contact problem and again we refer to
-Hintermueller, Ito, Kunisch: The primal-dual active set strategy as a semismooth newton method, SIAM J. OPTIM., 2003, Vol. 13, No. 3, pp. 865-888.
+Hintermueller, Ito, Kunisch: The primal-dual active set strategy as a semismooth Newton method, SIAM J. OPTIM., 2003, Vol. 13, No. 3, pp. 865-888.
Again the first case is for elastic and the second for plastic deformation.
\section{Formulation as a saddle point problem}
in step-41. The only difference consists in the fact that the contact area
is at the boundary instead of in the domain. But this has no further consequence
so that we refer to the documentation of step-41 with the only hint that
-$\mathcal{S}$ containts all the vertices at the contact boundary $\Gamma_C$ this
+$\mathcal{S}$ contains all the vertices at the contact boundary $\Gamma_C$ this
time.
\section{The primal-dual active set algorithm combined with the inexact semi smooth
\item[(0)] Initialize $\mathcal{A}_k$ and $\mathcal{F}_k$, such that
$\mathcal{S} = \mathcal{A}_k \cup \mathcal{F}_k$ and $\mathcal{A}_k \cap
\mathcal{F}_k = \emptyset$ and set $k = 1$. The start value $\hat U^0 :=
- P_{\mathcal{A}_k}(0)$ fullfills our obstacle condition.
+ P_{\mathcal{A}_k}(0)$ fulfills our obstacle condition.
\item[(1)] Assemble the Newton matrix $A := a'(\hat
U^{k-1};\varphi_p,\varphi_q)$ and the right-hand-side $F(\hat U^{k-1})$.
\item[(2)] Find the primal-dual pair $(\bar U^k,\Lambda^k)$ that satisfies
p\in\mathcal{A}.
\end{cases}$$\\
The matrix $B\in\mathbb{R}^{n\times m}$, $n>m$ describes the coupling of the
-bases for the displacements and lagrange multiplier (contact forces)
+bases for the displacements and Lagrange multiplier (contact forces)
and it is not quadratic in our situation since $\Lambda^k$ is only defined on
$\Gamma_C$. Due to the ansatz functions $\psi_i$ (scalar valued) of the
-lagrange multiplier are fullfilling the following biorthogonal condition (see Hüeber, Wohlmuth: A primal–dual active
+Lagrange multiplier are fulfilling the following biorthogonal condition (see Hüeber, Wohlmuth: A primal–dual active
set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg.
194, 2005, pp. 3147-3166)
$$ \int\limits_{\Gamma_C}\psi_i(x)\varphi_j(x)dx =
p,q\in\mathcal{S}\\
0I_3, & \text{if}\quad p\neq q,\quad p\textrm{ or }q\notin\mathcal{S}.
\end{cases}$$
-Here $I_3$ denotes the threedimensional identity matrix.
-In our programm we use the structure of a quadratic sparse for
+Here $I_3$ denotes the three-dimensional identity matrix.
+In our program we use the structure of a quadratic sparse for
$B\in\mathbb{R}^{n\times n}$ and for $\Lambda^k$ a vector with length $n$ where
$\Lambda^k_p = 0$ for $p\notin \mathcal{S}$.
The vector $G$ is defined by a suitable approximation $g_h$ of the gap $g$
g_{h,p}\int\limits_{\Gamma_C}\varphi_p(x)dx, & \text{if}\quad p\in\mathcal{S}\\
0, & \text{if}\quad p\notin\mathcal{S}.
\end{cases}$$\\
-Note that $G_p$ is a threedimensional vector and that again we applied the
-biorthogonal property of the lagrange multiplier ansatz functions to the
+Note that $G_p$ is a three-dimensional vector and that again we applied the
+biorthogonal property of the Lagrange multiplier ansatz functions to the
integral $\int\limits_{\Gamma_C}g_h(x)\varphi_p(x)dx$ with $g_h(x)=\sum\limits_i
g_{h,p}\varphi_p(x)$ (see the reference mentioned above).\\
Compared to step-41, step (1) is added but it should be clear from the sections
\section{Implementation}
-This tutorial is essentailly a mixture of step-40 and step-41 but instead of
+This tutorial is essentially a mixture of step-40 and step-41 but instead of
PETSc we let the Trilinos library deal with parallelizing the linear algebra
(like in step-32). Since we are trying to solve a similar problem like in
step-41 we will use the same methods but now in parallel.
Another difficulty is the handling of the different constraints from
-(the dirichlet conditons), the hanging nodes and the inequality condition that
+(the Dirichlet conditions), the hanging nodes and the inequality condition that
arises from the contact. For this purpose we create three objects of type
ConstraintMatrix.
Beside the ConstitutiveLaw class there is another new class called Input. This
class allows us to read in an obstacle from a file. In our example the file
-'obstacle\_file.dat' containts data which describe an Chinese, Japanese or
+'obstacle\_file.dat' contains data which describe an Chinese, Japanese or
Korean symbol for force or power. (See www.orientaloutpost.com/: ``This word can be used for motivation - it
can also mean power / motion / propulsion / force. It can be anything
internal or external that keeps you going. This is the safest way to express
To consider it the other way around if the deviator stress $\sigma^D$ is in a
norm bigger than the yield stress then $\sigma^D$ has to be projected back to the yield surface and there are plastic deformations which means $\lambda$
-would be positiv for that particular point. We refer that the stresses are
-computed by Hooke's law for isotorpic materials. You can find the description at the end of section 3. Else if the norm of the deviator stress tensor is smaller or equal the yield stress then $\lambda$ is zero and there are no plastic deformations in
+would be positive for that particular point. We refer that the stresses are
+computed by Hooke's law for isotropic materials. You can find the description at the end of section 3. Else if the norm of the deviator stress tensor is smaller or equal the yield stress then $\lambda$ is zero and there are no plastic deformations in
that point.
There the index $D$ denotes the deviator part of for example the stress where
The boundary of $\Omega$ separates as follows $\Gamma=\Gamma_D\bigcup\Gamma_C$ and $\Gamma_D\bigcap\Gamma_C=\emptyset$.
At the boundary $\Gamma_D$ we have zero Dirichlet conditions. $\Gamma_C$ denotes the potential contact boundary.
-The last two lines decribe the so-called Signorini contact conditions. If there is no contact the normal stress
+The last two lines describe the so-called Signorini contact conditions. If there is no contact the normal stress
@f{gather*} \sigma_n = \sigma n\cdot n@f}
is zero with the outward normal $n$. If there is contact ($u_n = g$) the tangential stress $\sigma_t = \sigma\cdot n - \sigma_n n$
vanishes, because we consider a frictionless situation and the normal stress is
L^2(\Omega,\mathbb{R}^{\textrm{dim}\times\textrm{dim}}_{\textrm{sym}}):\textrm{div}(\tau)\in L^2(\Omega,\mathbb{R}^{\textrm{dim}})\rbrace@f} and
@f{gather*}\Pi \Sigma:=\lbrace \tau\in \Sigma, \mathcal{F}(\tau)\leq 0\rbrace@f}
as the set of admissible stresses which is defined
-by a continious, convex flow function $\mathcal{F}$.
+by a continuous, convex flow function $\mathcal{F}$.
With the goal of deriving the dual formulation of the minimisation
-problem, we define a lagrange function:
+problem, we define a Lagrange function:
@f{gather*}L(\tau,\varphi) := E(\tau) + (\varphi, \textrm{div}(\tau)),\quad \lbrace\tau,\varphi\rbrace\in\Pi W^{\textrm{div}}\times V^+@f}
with
@f{gather*}V^+ := \lbrace u\in V: u_n\leq g \text{ on } \Gamma_C \rbrace@f}
@f{gather*}\tau^D := C\varepsilon^D(u^i).@f}
Remark that $a(.;.)$ is not differentiable in the common sense but it is
slantly differentiable like the function for the contact problem and again we refer to
-Hintermueller, Ito, Kunisch: The primal-dual active set strategy as a semismooth newton method, SIAM J. OPTIM., 2003, Vol. 13, No. 3, pp. 865-888.
+Hintermueller, Ito, Kunisch: The primal-dual active set strategy as a semismooth Newton method, SIAM J. OPTIM., 2003, Vol. 13, No. 3, pp. 865-888.
Again the first case is for elastic and the second for plastic deformation.
in step-41. The only difference consists in the fact that the contact area
is at the boundary instead of in the domain. But this has no further consequence
so that we refer to the documentation of step-41 with the only hint that
-$\mathcal{S}$ containts all the vertices at the contact boundary $\Gamma_C$ this
+$\mathcal{S}$ contains all the vertices at the contact boundary $\Gamma_C$ this
time.
<li> Initialize $\mathcal{A}_k$ and $\mathcal{F}_k$, such that
$\mathcal{S} = \mathcal{A}_k \cup \mathcal{F}_k$ and $\mathcal{A}_k \cap
\mathcal{F}_k = \emptyset$ and set $k = 1$. The start value $\hat U^0 :=
- P_{\mathcal{A}_k}(0)$ fullfills our obstacle condition.
+ P_{\mathcal{A}_k}(0)$ fulfills our obstacle condition.
<li> Assemble the Newton matrix $A := a'(\hat
U^{k-1};\varphi_p,\varphi_q)$ and the right-hand-side $F(\hat U^{k-1})$.
<li> Find the primal-dual pair $(\bar U^k,\Lambda^k)$ that satisfies
p\in\mathcal{A}.
\end{cases}@f}\\
The matrix $B\in\mathbb{R}^{n\times m}$, $n>m$ describes the coupling of the
-bases for the displacements and lagrange multiplier (contact forces)
+bases for the displacements and Lagrange multiplier (contact forces)
and it is not quadratic in our situation since $\Lambda^k$ is only defined on
$\Gamma_C$. Due to the ansatz functions $\psi_i$ (scalar valued) of the
-lagrange multiplier are fullfilling the following biorthogonal condition (see Hüeber, Wohlmuth: A primal–dual active
+Lagrange multiplier are fulfilling the following biorthogonal condition (see Hüeber, Wohlmuth: A primal–dual active
set strategy for non-linear multibody contact problems, Comput. Methods Appl. Mech. Engrg.
194, 2005, pp. 3147-3166)
@f{gather} \int\limits_{\Gamma_C}\psi_i(x)\varphi_j(x)dx =
p,q\in\mathcal{S}\\
0I_3, & \text{if}\quad p\neq q,\quad p\textrm{ or }q\notin\mathcal{S}.
\end{cases}@f}
-Here $I_3$ denotes the threedimensional identity matrix.
-In our programm we use the structure of a quadratic sparse for
+Here $I_3$ denotes the three-dimensional identity matrix.
+In our program we use the structure of a quadratic sparse for
$B\in\mathbb{R}^{n\times n}$ and for $\Lambda^k$ a vector with length $n$ where
$\Lambda^k_p = 0$ for $p\notin \mathcal{S}$.
The vector $G$ is defined by a suitable approximation $g_h$ of the gap $g$
0, & \text{if}\quad p\notin\mathcal{S}.
\end{cases}@f}
-Note that $G_p$ is a threedimensional vector and that again we applied the
-biorthogonal property of the lagrange multiplier ansatz functions to the
+Note that $G_p$ is a three-dimensional vector and that again we applied the
+biorthogonal property of the Lagrange multiplier ansatz functions to the
integral $\int\limits_{\Gamma_C}g_h(x)\varphi_p(x)dx$ with $g_h(x)=\sum\limits_i
g_{h,p}\varphi_p(x)$ (see the reference mentioned above).
step-41 we will use the same methods but now in parallel.
Another difficulty is the handling of the different constraints from
-(the dirichlet conditons), the hanging nodes and the inequality condition that
+(the Dirichlet conditions), the hanging nodes and the inequality condition that
arises from the contact. For this purpose we create three objects of type
ConstraintMatrix.
Beside the ConstitutiveLaw class there is another new class called Input. This
class allows us to read in an obstacle from a file. In our example the file
-'obstacle_file.dat' containts data which describe an Chinese, Japanese or
+'obstacle_file.dat' contains data which describe an Chinese, Japanese or
Korean symbol for force or power. (See http://www.orientaloutpost.com/ :
"This word can be used for motivation - it
can also mean power/motion/propulsion/force. It can be anything
For each adaptive refinement cycle the iterations end if the active set doesn't
change anymore and if the residual is accurate enough. In the tables
-at the end of each cycle you find informations about computing time and the
+at the end of each cycle you find information about computing time and the
number of calls of different parts of the program like Assembling or calculating
the residual. In the second cycle there are 12 calls for the residual function
and only 10 calls for Solving or Assembling what means that two damping steps were
Chinese symbol as obstacle it turns out that an elastic start solution for each
cycle results in fewer Newton iterations.
-In every refinement step you can observe that the acitve set - the contact points -
+In every refinement step you can observe that the active set - the contact points -
are iterated out at first. After that the Newton method has only to struggle with the
plasticity. For the finer meshes there is quadratic convergence not until the
last 4 or 5 Newton iterations.
</table>
The picture shows the adaptive refinement and as well how much a cell is
-plastified druing the contact with the ball. Remember that we consider the
+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. In the middle of the top -
where the mesh is finest - you can see the hollow caused by the ball. The blue
<h1>Possibilities for extensions</h1>
Extend the program from a static to a quasi-static problem, perhaps by choosing a
-backward-euler-scheme for the time discretization (for theoretical results see Frohne: FEM-Simulation
+backward-Euler-scheme for the time discretization (for theoretical results see Frohne: FEM-Simulation
der Umformtechnik metallischer Oberflächen im Mikrokosmos, Ph.D. thesis,
University of Siegen, Germany, 2011).
locally_relevant_dofs);
}
- // setup hanging nodes and dirichlet constraints
+ // setup hanging nodes and Dirichlet constraints
{
TimerOutput::Scope t(computing_timer, "Setup: constraints");
constraints_hanging_nodes.reinit(locally_relevant_dofs);
// This function defines the new ConstraintMatrix
// constraints_dirichlet_hanging_nodes. It contains
-// the dirichlet boundary values as well as the
+// the Dirichlet boundary values as well as the
// hanging nodes constraints.
template <int dim>
void
<li>
JL Guermond and R Pasquetti.
-<br> Entropy-based nonlinear viscosity for fourier approximations of
+<br> Entropy-based nonlinear viscosity for Fourier approximations of
conservation laws.
<br> <i>Comptes Rendus Mathematique</i>, 346(13-14):801-806, 2008.
// Unlike step-31, this step uses one more ConstraintMatrix object called
// darcy_preconditioner_constraints. This constraint object is used only for
// assembling the matrix for the Darcy preconditioner and includes hanging
- // node constrants as well as Dirichlet boundary value constraints for the
+ // node constraints as well as Dirichlet boundary value constraints for the
// pressure variable. We need this because we are building a Laplace matrix
// for the pressure as an approximation of the Schur complement) which is
// only positive definite if boundary conditions are applied.
// the run() function, the central one in this program.
//
// At the beginning of the function, we ask whether to solve the
- // pressure-velocity part by evaluating the posteriori criterion (see the
+ // pressure-velocity part by evaluating the a posteriori criterion (see the
// following function). If necessary, we will solve the pressure-velocity
// part using the GMRES solver with the Schur complement block
// preconditioner as is described in the introduction.
//
// With the exception of the startup code that loops back to the beginning
// of the function through the <code>goto start_time_iteration</code> label,
- // everything should be relatively straightforward. In any case, it mimicks
+ // everything should be relatively straightforward. In any case, it mimics
// the corresponding function in step-31.
template <int dim>
void TwoPhaseFlowProblem<dim>::run ()
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.
-One such straegy is the three-field formulation.
+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.
Additionally, the three-field formulation employed is valid for quasi-incompressible as well as compressible materials.
@f}
There are no derivatives of the pressure and dilatation (primary) variables present in the formulation.
-Thus the discontinuous finite element interpolation of the pressure ana dilatation yields a block
+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
discussed and demonstrated in Step-18 and others.
With over 93% of the time being spent in the linear solver, it is obvious that it may be necessary
to invest in a better solver for large three-dimensional problems.
-The SSOR preconditioner is not multi-threaded but is effective for this class of solid problems.
+The SSOR preconditioner is not multithreaded but is effective for this class of solid problems.
It may be beneficial to investigate the use of another solver such as those available through the Trilinos library.
}
protected:
- // Define constitutive model paramaters $\kappa$ (bulk modulus) and the
+ // Define constitutive model parameters $\kappa$ (bulk modulus) and the
// neo-Hookean model parameter $c_1$:
const double kappa;
const double c_1;
// \textrm{Grad}\ \mathbf{u}$: Since $I$ has data type SymmetricTensor,
// just writing <code>I + Grad_u_n</code> would convert the second
// argument to a symmetric tensor, perform the sum, and then cast the
- // result to a Tensor (i.e., the type of a possibly non-symmetric
+ // result to a Tensor (i.e., the type of a possibly nonsymmetric
// tensor). However, since <code>Grad_u_n</code> is nonsymmetric in
// general, the conversion to SymmetricTensor will fail. We can avoid this
// back and forth by converting $I$ to Tensor first, and then performing
- // the addition as between non-symmetric tensors:
+ // the addition as between nonsymmetric tensors:
void update_values (const Tensor<2, dim> &Grad_u_n,
const double p_tilde,
const double J_tilde)
determine_component_extractors();
// Several functions to assemble the system and right hand side matrices
- // using multi-threading. Each of them comes as a wrapper function, one
+ // using multithreading. Each of them comes as a wrapper function, one
// that is executed to do the work in the WorkStream model on one cell,
// and one that copies the work done on this one cell into the global
// object that represents it:
std::vector<PointHistory<dim> > quadrature_point_history;
// A description of the finite-element system including the displacement
- // polynomial degree, the degree-of-freedom handler, number of dof's per
+ // polynomial degree, the degree-of-freedom handler, number of DoFs per
// cell and the extractor objects used to retrieve information from the
// solution vectors:
const unsigned int degree;
block_component[J_component] = J_dof; // Dilatation
// The DOF handler is then initialised and we renumber the grid in an
- // efficient manner. We also record the number of DOF's per block.
+ // efficient manner. We also record the number of DOFs per block.
dof_handler_ref.distribute_dofs(fe);
DoFRenumbering::Cuthill_McKee(dof_handler_ref);
DoFRenumbering::component_wise(dof_handler_ref, block_component);
// @sect4{Solid::update_qph_incremental}
// As the update of QP information occurs frequently and involves a number of
-// expensive operations, we define a multi-threaded approach to distributing
+// expensive operations, we define a multithreaded approach to distributing
// the task across a number of CPU cores.
//
// To start this, we first we need to obtain the total solution as it stands
// @sect4{Solid::solve_linear_system}
// Solving the entire block system is a bit problematic as there are no
// contributions to the $\mathsf{\mathbf{k}}_{ \widetilde{J} \widetilde{J}}$
-// block, rendering it non-invertible.
+// block, rendering it noninvertible.
// Since the pressure and dilatation variables DOFs are discontinuous, we can
// condense them out to form a smaller displacement-only system which
// we will then solve and subsequently post-process to retrieve the
- Implementing the bilinear form, and in particular dealing with the
interface term, both in the matrix and the sparsity pattern.
- Implementing Dirichlet boundary conditions on the external and
- internal parts of the boundaryies
+ internal parts of the boundaries
$\partial\Omega_f,\partial\Omega_s$.
// The following classes do as their names suggest. The boundary values for
// the velocity are $\mathbf u=(0, \sin(\pi x))^T$ in 2d and $\mathbf u=(0,
// 0, \sin(\pi x)\sin(\pi y))^T$ in 3d, respectively. The remaining boundary
- // conditions for this problem are all homogenous and have been discussed in
+ // conditions for this problem are all homogeneous and have been discussed in
// the introduction. The right hand side forcing term is zero for both the
// fluid and the solid.
template <int dim>
// this end, we first have to set the active FE indices with the function
// immediately above, then distribute degrees of freedom, and then determine
// constraints on the linear system. The latter includes hanging node
- // constraints as usual, but also the inhomogenous boundary values at the
+ // constraints as usual, but also the inhomogeneous boundary values at the
// top fluid boundary, and zero boundary values along the perimeter of the
// solid subdomain.
template <int dim>
// To integrate the enriched elements we have to find the geometrical
// decomposition of the original element in subelements. The subelements are
// used to integrate the elements on both sides of the discontinuity. The
-// disontinuity line is approximated by a piece-wise linear interpolation
+// discontinuity line is approximated by a piece-wise linear interpolation
// between the intersection of the discontinuity with the edges of the
// elements. The vector level_set_values has the values of the level set
// function at the vertices of the elements. From these values can be found by
// the sign of the level set function at the 4 nodes of the elements can
// be positive + or negative - depending on the sign of the level set
- // function we have the folloing three classes of decomposition type 1:
+ // function we have the following three classes of decomposition type 1:
// ++++, ---- type 2: -+++, +-++, ++-+, +++-, +---, -+--, --+-, ---+ type
// 3: +--+, ++--, +-+-, -++-, --++, -+-+
indices of the DoFs that it is constrained to. Then, we have another
variable <code>constraint_indicator</code> at hand that holds, for
each cell, the local indices of DoFs that are constrained as well as
-the identifer of the type of constraint. Actually, you will not see
+the identifier of the type of constraint. Actually, you will not see
these data structures in the example program since the class
<code>FEEvaluationGL</code> takes care of the constraints without user
interaction.
<h3>Parallel run in 3D</h3>
-To demonstrate how the example scales for a parallel run and to demsonstrate
+To demonstrate how the example scales for a parallel run and to demonstrate
that hanging node constraints can be handled in an efficient way, we run the
example in 3D with $\mathcal{Q}_4$ elements. First, we run it on a notebook
with 2 cores (Sandy Bridge CPU) at 2.7 GHz.
@endcode
We observe a considerable speedup over the notebook (16 cores versus 2 cores;
-nonethess, one notebook core is considerably faster than one core of the
+nonetheless, one notebook core is considerably faster than one core of the
cluster because of a newer processor architecture). If we run the same program
on 4 nodes with 8 threads on each node, we get:
@code
// dimensions of extent $[-15,15]$. We refine the mesh more in the center of
// the domain since the solution is concentrated there. We first refine all
// cells whose center is within a radius of 11, and then refine once more
- // for a radius 6. This is simple ad-hoc refinement could be done better by
+ // for a radius 6. This simple ad hoc refinement could be done better by
// adapting the mesh to the solution using error estimators during the time
// stepping as done in other example programs, and using
// parallel::distributed::SolutionTransfer to transfer the solution to the
<h3>Modifying a Mesh</h3>
-After aquiring one (or several) meshes in the ways described above, there are
+After acquiring one (or several) meshes in the ways described above, there are
many ways to manipulate them before using them in a finite element
computation.
<h1>Results</h1>
The program produces a series of <code>.eps</code> files of the
-Triangulations. The methods are discussed above.
+triangulations. The methods are discussed above.
<h3>Next steps: Curved boundaries</h3>
*/
-// parallel geometric multi-grid. work in progress!
+// parallel geometric multigrid. work in progress!
// As discussed in the introduction, most of
// this program is copied almost verbatim
#include <deal.II/lac/trilinos_precondition.h>
// These, now, are the include necessary for
-// the multi-level methods. The first two
+// the multilevel methods. The first two
// declare classes that allow us to enumerate
// degrees of freedom not only on the finest
// mesh level, but also on intermediate
// Now for the things that concern the
// multigrid data structures. First, we
- // resize the multi-level objects to hold
+ // resize the multilevel objects to hold
// matrices and sparsity patterns for every
// level. The coarse level is zero (this is
// mandatory right now but may change in a
// @sect4{LaplaceProblem::assemble_system}
// The following function assembles the
- // linear system on the finesh level of the
+ // linear system on the finest level of the
// mesh. It is almost exactly the same as in
// step-6, with the exception that we don't
// eliminate hanging nodes and boundary
// over all existing cells instead of just
// the active ones, and the results must be
// entered into the correct matrix. Note also
- // that since we only do multi-level
+ // that since we only do multilevel
// preconditioning, no right-hand side needs
// to be assembled here.
//
// level that are located on interfaces between adaptively refined
// levels, and those that lie on the interface but also on the
// exterior boundary of the domain. As in many other parts of the
- // library, we do this by using boolean masks, i.e. vectors of
- // booleans each element of which indicates whether the
+ // library, we do this by using Boolean masks, i.e. vectors of
+ // Booleans each element of which indicates whether the
// corresponding degree of freedom index is an interface DoF or
// not. The <code>MGConstraints</code> already computed the
// information for us when we called initialize in
MGTransferPrebuilt<vector_t> mg_transfer(hanging_node_constraints, mg_constrained_dofs);
// Now the prolongation matrix has to be built. This matrix needs
// to take the boundary values on each level into account and
- // needs to know about the indices at the refinement egdes. The
+ // needs to know about the indices at the refinement edges. The
// <code>MGConstraints</code> knows about that so pass it as an
// argument.
mg_transfer.build_matrices(mg_dof_handler);
// difference is the use of 3 different sets of <code>DoFHandler</code> and FE
// objects, along with the <code>ChunkSparseMatrix</code> and the
// corresponding solutions vectors. We also use WorkStream to enable a
-// multi-threaded local solution process which exploits the embarrassingly
+// multithreaded local solution process which exploits the embarrassingly
// parallel nature of the local solver. For WorkStream, we define the local
// operations on a cell and a copy function into the global matrix and
// vector. We do this once for the assembly (which is run twice, once when we
Vector<double> system_rhs;
// As stated in the introduction, HDG solutions can be post-processed to
- // attain superconvegence rates of $\mathcal{O}(h^{p+2})$. The
+ // attain superconvergence rates of $\mathcal{O}(h^{p+2})$. The
// post-processed solution is a discontinuous finite element solution
// representing the primal variable on the interior of each cell. We define
// a FE type of degree $p+1$ to represent this post-processed solution,
#include <deal.II/grid/grid_refinement.h>
// Finally, we need a simple way to actually compute the refinement indicators
-// based on some error estimat. While in general, adaptivity is very
+// based on some error estimate. While in general, adaptivity is very
// problem-specific, the error indicator in the following file often yields
// quite nicely adapted grids for a wide class of problems.
#include <deal.II/numerics/error_estimator.h>
// After setting up all the degrees of freedoms, here are now the
// differences compared to step-5, all of which are related to constraints
- // associated with the hanging nodes. In the class desclaration, we have
+ // associated with the hanging nodes. In the class declaration, we have
// already allocated space for an object <code>constraints</code> that will
// hold a list of these constraints (they form a matrix, which is reflected
// in the name of the class, but that is immaterial for the moment). Now we
// approximately the same string as would be generated if the exception was
// thrown using the <code>Assert</code> macro. You have seen the output of
// such an exception in the previous example, and you then know that it
- // contains the file and line number of where the exception occured, and
+ // contains the file and line number of where the exception occurred, and
// some other information. This is also what the following statements would
// print.
//
that can be obtained by assuming that the body has infinite extent in
$z$-direction and that the solution function does not depend on the $z$
coordinate. On the other hand, there are equations for two-dimensional models
-of elasticity; see for example the
wikipedia article on <a
+of elasticity; see for example the
Wikipedia article on <a
href="http://en.wikipedia.org/wiki/Infinitesimal_strain_theory#Special_cases">plane
strain</a>, <a
href="http://en.wikipedia.org/wiki/Antiplane_shear">antiplane shear</a> and <a
// @sect4{ElasticProblem::setup_system}
- // Setting up the system of equations is identitical to the function used in
+ // Setting up the system of equations is identical to the function used in
// the step-6 example. The <code>DoFHandler</code> class and all other
// classes used here are fully aware that the finite element we want to use
// is vector-valued, and take care of the vector-valuedness of the finite
#include <deal.II/fe/fe_q.h>
#include <deal.II/grid/grid_out.h>
-// The following two files provide classes and information for multi-threaded
+// The following two files provide classes and information for multithreaded
// programs. In the first one, the classes and functions are declared which we
// need to start new threads and to wait for threads to return (i.e. the
// <code>Thread</code> class and the <code>new_thread</code> functions). The
// information, we could use the value of the global variable
// <code>multithread_info.n_cpus</code>, which is determined at start-up
// time of your program automatically. (Note that if the library was not
- // configured for multi-threading, then the number of CPUs is set to one.)
+ // configured for multithreading, then the number of CPUs is set to one.)
// However, sometimes there might be reasons to use another value. For
// example, you might want to use less processors than there are in your
// system in order not to use too many computational resources. On the
// can be omitted. (However, you still need to write the angle brackets,
// even if they are empty.)
//
- // If you did not configure for multi-threading, then the
+ // If you did not configure for multithreading, then the
// <code>new_thread</code> function that is supposed to start a new thread
// in parallel only executes the function which should be run in parallel,
// waits for it to return (i.e. the function is executed sequentially),
// the same size. Each thread will then assemble the local contributions
// of the cells within its chunk and transfer these contributions to the
// global matrix. As splitting a range of cells is a rather common task
- // when using multi-threading, there is a function in the
+ // when using multithreading, there is a function in the
// <code>Threads</code> namespace that does exactly this. In fact, it does
// this not only for a range of cell iterators, but for iterators in
// general, so you could use it for <code>std::vector::iterator</code> or
n_threads);
// Finally, for each of the chunks of iterators we have computed, start
- // one thread (or if not in multi-thread mode: execute assembly on these
+ // one thread (or if not in multithread mode: execute assembly on these
// chunks sequentially). This is done using the following sequence of
// function calls:
for (unsigned int thread=0; thread<n_threads; ++thread)
thread_ranges[thread].first,
thread_ranges[thread].second);
// The reasons and internal workings of these functions can be found in
- // the report on the subject of multi-threading, which is available online
+ // the report on the subject of multithreading, which is available online
// as well. Suffice it to say that we create a new thread that calls the
// <code>assemble_system_interval</code> function on the present object
// (the <code>this</code> pointer), with the arguments following in the
// container, which just calls <code>join</code> on each of the thread
// objects it stores.
//
- // Again, if the library was not configured to use multi-threading, then
+ // Again, if the library was not configured to use multithreading, then
// no threads can run in parallel and the function returns immediately.
threads.join_all ();
// i.e. other threads can now enter into the protected section by
// acquiring the lock. Two final notes are in place here, however:
//
- // 1. If the library was not configured for multi-threading, then
+ // 1. If the library was not configured for multithreading, then
// there can't be parallel threads and there is no need to
// synchronize. Thus, the <code>lock</code> and <code>release</code>
// functions are no-ops, i.e. they return without doing anything.
// can determine an approximation of the gradient for the present
// cell, then we need to have passed over vectors <code>y</code> which
// span the whole space, otherwise we would not have all components of
- // the gradient. This is indicated by the invertability of the matrix.
+ // the gradient. This is indicated by the invertibility of the matrix.
//
// If the matrix should not be invertible, this means that the present
// cell had an insufficient number of active neighbors. In contrast to
#cmakedefine DEAL_II_WITH_P4EST
#ifdef DEAL_II_WITH_P4EST
# define DEAL_II_USE_P4EST
+
+# define DEAL_II_P4EST_VERSION_MAJOR @P4EST_VERSION_MAJOR@
+# define DEAL_II_P4EST_VERSION_MINOR @P4EST_VERSION_MINOR@
+# define DEAL_II_P4EST_VERSION_SUBMINOR @P4EST_VERSION_SUBMINOR@
+# define DEAL_II_P4EST_VERSION_PATCH @P4EST_VERSION_PATCH@
+
+# define DEAL_II_P4EST_VERSION_GTE(major,minor,subminor,patch) \
+ ((DEAL_II_P4EST_VERSION_MAJOR * 1000000 + \
+ DEAL_II_P4EST_VERSION_MINOR * 10000 + \
+ DEAL_II_P4EST_VERSION_SUBMINOR * 100 + \
+ DEAL_II_P4EST_VERSION_PATCH) \
+ >= \
+ (major)*1000000 + (minor)*10000 + (subminor)*100 + (patch))
+#else
+ // p4est up to 0.3.4.1 didn't define P4EST_VERSION_*. since
+ // we didn't supports anything before 0.3.4, we assume 0.3.4
+ // This means that we can't use the new features in 0.3.4.1
+# define DEAL_II_P4EST_VERSION_GTE(major,minor,subminor,patch) \
+ ((0 * 1000000 + \
+ 3 * 10000 + \
+ 4 * 100 + \
+ 0) \
+ >= \
+ (major)*1000000 + (minor)*10000 + (subminor)*100 + (patch))
+
#endif
* these variables at a
* couple places anyway.
*/
- std::vector<types::global_dof_index> coarse_cell_to_p4est_tree_permutation;
- std::vector<types::global_dof_index> p4est_tree_to_coarse_cell_permutation;
+ std::vector<unsigned int> coarse_cell_to_p4est_tree_permutation;
+ std::vector<unsigned int> p4est_tree_to_coarse_cell_permutation;
/**
* dummy settings
const unsigned int local_index,
const dealii::internal::int2type<1> &)
{
- return dof_handler.levels[obj_level]->dof_object.
- get_dof_index (dof_handler,
- obj_index,
+ return dof_handler.levels[obj_level]->
+ get_dof_index (obj_index,
fe_index,
- local_index,
- obj_level);
+ local_index);
}
const dealii::internal::int2type<1> &,
const types::global_dof_index global_index)
{
- dof_handler.levels[obj_level]->dof_object.
- set_dof_index (dof_handler,
- obj_index,
+ dof_handler.levels[obj_level]->
+ set_dof_index (obj_index,
fe_index,
local_index,
- global_index,
- obj_level);
+ global_index);
}
{
return dof_handler.faces->lines.
get_dof_index (dof_handler,
- obj_index,
+ obj_index,
fe_index,
local_index,
obj_level);
const unsigned int local_index,
const dealii::internal::int2type<2> &)
{
- return dof_handler.levels[obj_level]->dof_object.
- get_dof_index (dof_handler,
- obj_index,
+ return dof_handler.levels[obj_level]->
+ get_dof_index (obj_index,
fe_index,
- local_index,
- obj_level);
+ local_index);
}
const dealii::internal::int2type<2> &,
const types::global_dof_index global_index)
{
- dof_handler.levels[obj_level]->dof_object.
- set_dof_index (dof_handler,
- obj_index,
+ dof_handler.levels[obj_level]->
+ set_dof_index (obj_index,
fe_index,
local_index,
- global_index,
- obj_level);
+ global_index);
}
const unsigned int local_index,
const dealii::internal::int2type<3> &)
{
- return dof_handler.levels[obj_level]->dof_object.
- get_dof_index (dof_handler,
- obj_index,
+ return dof_handler.levels[obj_level]->
+ get_dof_index (obj_index,
fe_index,
- local_index,
- obj_level);
+ local_index);
}
const dealii::internal::int2type<3> &,
const types::global_dof_index global_index)
{
- dof_handler.levels[obj_level]->dof_object.
- set_dof_index (dof_handler,
- obj_index,
+ dof_handler.levels[obj_level]->
+ set_dof_index (obj_index,
fe_index,
local_index,
- global_index,
- obj_level);
+ global_index);
}
const unsigned int fe_index,
const dealii::internal::int2type<1> &)
{
- return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
+ return dof_handler.levels[obj_level]->fe_index_is_active(obj_index,
+ fe_index);
}
template <int spacedim>
static
unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
+ n_active_fe_indices (const dealii::hp::DoFHandler<1,spacedim> &,
const unsigned int obj_level,
const unsigned int obj_index,
const dealii::internal::int2type<1> &)
{
- return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
- obj_index);
+ // on a cell, the number of active elements is one
+ return 1;
}
const unsigned int n,
const dealii::internal::int2type<1> &)
{
- return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
+ Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
+ return dof_handler.levels[obj_level]->active_fe_index (obj_index);
}
const unsigned int fe_index,
const dealii::internal::int2type<2> &)
{
- return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
+ return dof_handler.levels[obj_level]->fe_index_is_active(obj_index,
+ fe_index);
}
template <int spacedim>
static
unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
+ n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &,
const unsigned int obj_level,
const unsigned int obj_index,
const dealii::internal::int2type<2> &)
{
- return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
- obj_index);
+ // on a cell, the number of active elements is one
+ return 1;
}
const unsigned int n,
const dealii::internal::int2type<2> &)
{
- return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
+ Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
+ return dof_handler.levels[obj_level]->active_fe_index (obj_index);
}
const unsigned int fe_index,
const dealii::internal::int2type<3> &)
{
- return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
- obj_index,
- fe_index,
- obj_level);
+ return dof_handler.levels[obj_level]->fe_index_is_active(obj_index,
+ fe_index);
}
template <int spacedim>
static
unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
+ n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &,
const unsigned int obj_level,
const unsigned int obj_index,
const dealii::internal::int2type<3> &)
{
- return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
- obj_index);
+ // on a cell, the number of active elements is one
+ return 1;
}
const unsigned int n,
const dealii::internal::int2type<3> &)
{
- return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
- obj_level,
- obj_index,
- n);
+ Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
+ return dof_handler.levels[obj_level]->active_fe_index (obj_index);
}
/**
types::global_dof_index *cache = &accessor.dof_handler->levels[accessor.level()]
->cell_dof_indices_cache[accessor.present_index *
accessor.get_fe().dofs_per_cell];
- for ( ; local_values_begin != local_values_end; ++local_values_begin, ++cache)
- *local_values_begin = values(*cache);
+
+ values.extract_subvector_to (cache,
+ cache + accessor.get_fe().dofs_per_cell,
+ local_values_begin);
}
/**
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
get_dof_indices (accessor, local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i, ++local_values_begin)
- *local_values_begin = values(local_dof_indices[i]);
+ values.extract_subvector_to (local_dof_indices.begin(),
+ local_dof_indices.end(),
+ local_values_begin);
}
{
Assert (static_cast<unsigned int>(accessor.level()) < accessor.dof_handler->levels.size(),
ExcMessage ("DoFHandler not initialized"));
- Assert (static_cast<std::vector<unsigned int>::size_type>(accessor.present_index) <
- accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size (),
- ExcIndexRange (accessor.present_index, 0,
- accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size ()));
+
return accessor.dof_handler->levels[accessor.level()]
- ->active_fe_indices[accessor.present_index];
+ ->active_fe_index(accessor.present_index);
}
Assert (static_cast<unsigned int>(accessor.level()) <
accessor.dof_handler->levels.size(),
ExcMessage ("DoFHandler not initialized"));
- Assert (static_cast<std::vector<unsigned int>::size_type>(accessor.present_index) <
- accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size (),
- ExcIndexRange (accessor.present_index, 0,
- accessor.dof_handler->levels[accessor.level()]->active_fe_indices.size ()));
+
accessor.dof_handler->levels[accessor.level()]
- ->active_fe_indices[accessor.present_index] = i;
+ ->set_active_fe_index (accessor.present_index, i);
}
* this class should not be used for this
* context.
*
+ * This class allows to query the
+ * value, gradient and divergence of
+ * (components of) shape functions
+ * and solutions representing
+ * vectors. The
+ * gradient of a vector
+ * $d_{k}, 0\le k<\text{dim}$ is
+ * defined as
+ * $S_{ij} = \frac{\partial d_{i}}{\partial x_j},
+ * 0\le i,j<\text{dim}$.
+ *
* You get an object of this type if you
* apply a FEValuesExtractors::Vector to an
* FEValues, FEFaceValues or
* <code>dim</code> components of the
* finite element, the gradient is a
* <code>Tensor@<2,spacedim@></code>.
+ *
+ * See the general documentation of this class for how exactly
+ * the gradient of a vector is defined.
*/
typedef dealii::Tensor<2,spacedim> gradient_type;
* function and quadrature point
* selected by the arguments.
*
+ * See the general documentation of this class for how exactly
+ * the gradient of a vector is defined.
+ *
* @note The meaning of the arguments
* is as documented for the value()
* function.
*
* This function tries to match all faces belonging to the first
* boundary with faces belonging to the second boundary with the help
- * of orthogonal_equality.
+ * of orthogonal_equality().
+ *
+ * The bitset that is returned together with the second face encodes the
+ * _relative_ orientation of the first face with respect to the second
+ * face, see the documentation of orthogonal_equality for further details.
*
* The @p offset is a vector tangential to the faces that is added to the
* location of vertices of the 'first' boundary when attempting to match
const typename identity<FaceIterator>::type &end,
const types::boundary_id b_id1,
const types::boundary_id b_id2,
- int direction,
+ const int direction,
const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
/**
* Same function as above, but accepts a Triangulation or DoFHandler
- * object @p dof_handler instead of an explicit face iterator range.
+ * object @p container (a container is a collection of objects, here a
+ * collection of cells) instead of an explicit face iterator range.
*
* This function will collect periodic face pairs on the highest (i.e.
* coarsest) mesh level.
*/
template<typename DH>
std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
- collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/
+ collect_periodic_face_pairs (const DH &container,
const types::boundary_id b_id1,
const types::boundary_id b_id2,
- int direction,
+ const int direction,
const dealii::Tensor<1,DH::space_dimension> &offset);
std::map<typename DH::face_iterator, typename DH::face_iterator>
collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/
const types::boundary_id b_id,
- int direction,
+ const int direction,
const dealii::Tensor<1,DH::space_dimension> &offset);
#include <deal.II/base/config.h>
-#include <deal.II/hp/dof_objects.h>
+#include <deal.II/base/exceptions.h>
#include <vector>
+
DEAL_II_NAMESPACE_OPEN
+namespace hp
+{
+ template <int, int> class DoFHandler;
+}
+
+
namespace internal
{
namespace hp
{
- template <int dim>
- class DoFLevel;
+ namespace DoFHandler
+ {
+ struct Implementation;
+ }
+ }
+}
+
+namespace internal
+{
+ namespace hp
+ {
/**
- * Store the indices of the degrees of freedom that are located on
- * objects of dimension @p structdim.
- *
- * The things we store here are very similar to what is stored in the
- * internal::DoFHandler::DoFLevel class hierarchy (see there for more
- * information, in particular on the layout of the class hierarchy,
- * and the use of file names). There are two main
- * differences, discussed in the following subsections. In addition to
- * the data already stored by the internal::DoFHandler::DoFLevel
- * classes, we also have to store which finite element each cell
- * uses.
- *
- *
- * <h4>Offset computations</h4>
- *
- * For hp methods, not all cells may use the same finite element, and
- * it is consequently more complicated to determine where the DoF
- * indices for a given line, quad, or hex are stored. As described in
- * the documentation of the internal::DoFHandler::DoFLevel class, we
- * can compute the location of the first line DoF, for example, by
- * calculating the offset as <code>line_index *
- * dof_handler.get_fe().dofs_per_line</code>. This of course doesn't
- * work any more if different lines may have different numbers of
- * degrees of freedom associated with them. Consequently, rather than
- * using this simple multiplication, each of the lines.dofs, quads.dofs,
- * and hexes.dofs arrays has an associated array lines.dof_offsets,
- * quads.dof_offsets, and hexes.dof_offsets. The data corresponding to a
- * line then starts at index <code>line_dof_offsets[line_index]</code>
- * within the <code>line_dofs</code> array.
- *
- *
- * <h4>Multiple data sets per object</h4>
- *
- * If an object corresponds to a cell, the global dof indices of this
- * cell are stored at the location indicated above in sequential
- * order.
- *
- * However, if two adjacent cells use different finite elements, then
- * the face that they share needs to store DoF indices for both
- * involved finite elements. While faces therefore have to have at
- * most two sets of DoF indices, it is easy to see that vertices for
- * example can have as many sets of DoF indices associated with them
- * as there are adjacent cells, and the same holds for lines in 3d.
- *
- * Consequently, for objects that have a lower dimensionality than
- * cells, we have to store a map from the finite element index to the
- * set of DoF indices associated. Since real sets are typically very
- * inefficient to store, and since most of the time we expect the
- * number of individual keys to be small (frequently, adjacent cells
- * will have the same finite element, and only a single entry will
- * exist in the map), what we do is instead to store a linked list. In
- * this format, the first entry starting at position
- * <code>lines.dofs[lines.dof_offsets[line_index]]</code> will denote
- * the finite element index of the set of DoF indices following; after
- * this set, we will store the finite element index of the second set
- * followed by the corresponding DoF indices; and so on. Finally, when
- * all finite element indices adjacent to this object have been
- * covered, we write a -1 to indicate the end of the list.
- *
- * Access to this kind of data, as well as the distinction between
- * cells and objects of lower dimensionality are encoded in the
- * accessor functions, DoFObjects::set_dof_index() and
- * DoFLevel::get_dof_index() They are able to traverse this
- * list and pick out or set a DoF index given the finite element index
- * and its location within the set of DoFs corresponding to this
- * finite element.
- *
- *
- * @ingroup hp
- * @author Wolfgang Bangerth, 1998, 2006, Oliver Kayser-Herold 2003.
+ * This is the class that stores the degrees of freedom on cells in a hp
+ * hierarchy. Compared to faces and edges, the task here is simple since
+ * each cell can only have a single active finite element index. Consequently,
+ * all we need is one long array with DoF indices and one array of offsets
+ * where each cell's indices start within the array of indices. This is in
+ * contrast to the DoFObjects class where each face or edge may have more than
+ * one associated finite element with corresponding degrees of freedom.
*/
template <int dim>
class DoFLevel
{
- public:
+ private:
/**
- * Indices specifying the finite
- * element of hp::FECollection to use
- * for the different cells on the current level. The
- * meaning what a cell is, is
- * dimension specific, therefore also
- * the length of this vector depends
- * on the dimension: in one dimension,
- * the length of this vector equals
- * the length of the @p lines vector,
- * in two dimensions that of the @p
- * quads vector, etc. The vector stores one element per cell
- * since the actiev_fe_index is unique for cells.
+ * Indices specifying the finite element of hp::FECollection to
+ * use for the different cells on the current level. The vector
+ * stores one element per cell since the active_fe_index is
+ * unique for cells.
+ *
+ * If a cell is not active on the level corresponding to the
+ * current object (i.e., it has children on higher levels) then
+ * it does not have an associated fe index and we store
+ * an invalid fe index marker instead.
*/
std::vector<unsigned int> active_fe_indices;
/**
- * Store the dof-indices and related data of
- * the cells on the current level corresponding to this object.
+ * Store the start index for the degrees of freedom of each
+ * object in the @p dofs array. If the cell corresponding to
+ * a particular index in this array is not active on this level,
+ * then we do not store any DoFs for it. In that case, the offset
+ * we store here must be an invalid number and indeed we store
+ * <code>(std::vector<types::global_dof_index>::size_type)(-1)</code>
+ * for it.
+ *
+ * The type we store is then obviously the type the @p dofs array
+ * uses for indexing.
+ */
+ std::vector<std::vector<types::global_dof_index>::size_type> dof_offsets;
+
+ /**
+ * Store the global indices of the degrees of freedom.
+ * information. The dof_offsets field determines where each
+ * (active) cell's data is stored.
+ */
+ std::vector<types::global_dof_index> dofs;
+
+ public:
+
+ /**
+ * Set the global index of
+ * the @p local_index-th
+ * degree of freedom located
+ * on the object with number @p
+ * obj_index to the value
+ * given by @p global_index. The @p
+ * dof_handler argument is
+ * used to access the finite
+ * element that is to be used
+ * to compute the location
+ * where this data is stored.
+ *
+ * The third argument, @p
+ * fe_index, denotes which of
+ * the finite elements
+ * associated with this
+ * object we shall
+ * access. Refer to the
+ * general documentation of
+ * the internal::hp::DoFLevel
+ * class template for more
+ * information.
+ */
+ void
+ set_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const types::global_dof_index global_index);
+
+ /**
+ * Return the global index of
+ * the @p local_index-th
+ * degree of freedom located
+ * on the object with number @p
+ * obj_index. The @p
+ * dof_handler argument is
+ * used to access the finite
+ * element that is to be used
+ * to compute the location
+ * where this data is stored.
+ *
+ * The third argument, @p
+ * fe_index, denotes which of
+ * the finite elements
+ * associated with this
+ * object we shall
+ * access. Refer to the
+ * general documentation of
+ * the internal::hp::DoFLevel
+ * class template for more
+ * information.
+ */
+ types::global_dof_index
+ get_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index) const;
+
+ /**
+ * Return the fe_index of the
+ * active finite element
+ * on this object.
*/
- internal::hp::DoFObjects<dim> dof_object;
+ unsigned int
+ active_fe_index (const unsigned int obj_index) const;
+
+ /**
+ * Check whether a given
+ * finite element index is
+ * used on the present
+ * object or not.
+ */
+ bool
+ fe_index_is_active (const unsigned int obj_index,
+ const unsigned int fe_index) const;
+
+ /**
+ * Set the fe_index of the
+ * active finite element
+ * on this object.
+ */
+ void
+ set_active_fe_index (const unsigned int obj_index,
+ const unsigned int fe_index);
/**
* Determine an estimate for the
* of this object.
*/
std::size_t memory_consumption () const;
+
+ /**
+ * Make hp::DoFHandler and its auxiliary class a friend since it
+ * is the class that needs to create these data structures.
+ */
+ template <int, int> friend class dealii::hp::DoFHandler;
+ friend struct dealii::internal::hp::DoFHandler::Implementation;
};
+
+
+ // -------------------- template functions --------------------------------
+
+ template <int dim>
+ inline
+ types::global_dof_index
+ DoFLevel<dim>::
+ get_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index) const
+ {
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (fe_index == active_fe_indices[obj_index],
+ ExcMessage ("FE index does not match that of the present cell"));
+ return dofs[dof_offsets[obj_index]+local_index];
+ }
+
+
+
+ template <int dim>
+ inline
+ void
+ DoFLevel<dim>::
+ set_dof_index (const unsigned int obj_index,
+ const unsigned int fe_index,
+ const unsigned int local_index,
+ const types::global_dof_index global_index)
+ {
+ Assert (obj_index < dof_offsets.size(),
+ ExcIndexRange (obj_index, 0, dof_offsets.size()));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (fe_index == active_fe_indices[obj_index],
+ ExcMessage ("FE index does not match that of the present cell"));
+ dofs[dof_offsets[obj_index]+local_index] = global_index;
+ }
+
+
+
+ template <int dim>
+ inline
+ unsigned int
+ DoFLevel<dim>::
+ active_fe_index (const unsigned int obj_index) const
+ {
+ Assert (obj_index < active_fe_indices.size(),
+ ExcIndexRange (obj_index, 0, active_fe_indices.size()));
+
+ return active_fe_indices[obj_index];
+ }
+
+
+
+ template <int dim>
+ inline
+ bool
+ DoFLevel<dim>::
+ fe_index_is_active (const unsigned int obj_index,
+ const unsigned int fe_index) const
+ {
+ const unsigned int invalid_fe_index = numbers::invalid_unsigned_int;
+ Assert ((fe_index != invalid_fe_index),
+ ExcMessage ("You need to specify a FE index when working "
+ "with hp DoFHandlers"));
+
+ // make sure we are on an
+ // object for which DoFs have
+ // been allocated at all
+ Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ ExcMessage ("You are trying to access degree of freedom "
+ "information for an object on which no such "
+ "information is available"));
+
+ Assert (obj_index < active_fe_indices.size(),
+ ExcInternalError());
+ return (fe_index == active_fe_indices[obj_index]);
+ }
+
+
+ template <int dim>
+ inline
+ void
+ DoFLevel<dim>::
+ set_active_fe_index (const unsigned int obj_index,
+ const unsigned int fe_index)
+ {
+ Assert (obj_index < active_fe_indices.size(),
+ ExcIndexRange (obj_index, 0, active_fe_indices.size()));
+
+ active_fe_indices[obj_index] = fe_index;
+ }
+
} // namespace hp
} // namespace internal
/**
* Store the indices of the degrees of freedom which are located on
- * objects of dimension @p dim.
+ * objects of dimension @p structdim < dim, i.e., for faces or edges
+ * of cells. This is opposed to the internal::hp::DoFLevels class
+ * that stores the DoF indices on cells.
*
* The things we store here is very similar to what is stored in the
* internal::DoFHandler::DoFObjects classes (see there for more
*
* <h4>Multiple data sets per object</h4>
*
- * If an object corresponds to a cell, the global dof indices of this
- * cell are stored at the location indicated above in sequential
- * order.
- *
- * However, if two adjacent cells use different finite elements, then
+ * If two adjacent cells use different finite elements, then
* the face that they share needs to store DoF indices for both
* involved finite elements. While faces therefore have to have at
- * most two sets of DoF indices, it is easy to see that vertices for
- * example can have as many sets of DoF indices associated with them
- * as there are adjacent cells, and the same holds for lines in 3d.
+ * most two sets of DoF indices, it is easy to see that edges and
+ * vertices can have as many sets of DoF indices associated with them
+ * as there are adjacent cells.
*
* Consequently, for objects that have a lower dimensionality than
* cells, we have to store a map from the finite element index to the
* Access to this kind of data, as well as the distinction between
* cells and objects of lower dimensionality are encoded in the
* accessor functions, DoFObjects::set_dof_index() and
- * DoFLevel::get_dof_index() They are able to traverse this
+ * DoFLevel::get_dof_index(). They are able to traverse this
* list and pick out or set a DoF index given the finite element index
* and its location within the set of DoFs corresponding to this
* finite element.
* @ingroup hp
* @author Tobias Leicht, 2006
*/
- template <int dim>
+ template <int structdim>
class DoFObjects
{
public:
* class template for more
* information.
*/
- template <int dimm, int spacedim>
+ template <int dim, int spacedim>
void
- set_dof_index (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ set_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index,
const unsigned int fe_index,
const unsigned int local_index,
* class template for more
* information.
*/
- template <int dimm, int spacedim>
+ template <int dim, int spacedim>
types::global_dof_index
- get_dof_index (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ get_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index,
const unsigned int fe_index,
const unsigned int local_index,
* been distributed and zero
* is returned.
*/
- template <int dimm, int spacedim>
+ template <int dim, int spacedim>
unsigned int
- n_active_fe_indices (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ n_active_fe_indices (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index) const;
/**
* n-th active finite element
* on this object.
*/
- template <int dimm, int spacedim>
+ template <int dim, int spacedim>
types::global_dof_index
- nth_active_fe_index (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ nth_active_fe_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_level,
const unsigned int obj_index,
const unsigned int n) const;
* used on the present
* object or not.
*/
- template <int dimm, int spacedim>
+ template <int dim, int spacedim>
bool
- fe_index_is_active (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ fe_index_is_active (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index,
const unsigned int fe_index,
const unsigned int obj_level) const;
};
-// --------------------- inline and template functions ------------------
+ // --------------------- inline and template functions ------------------
- template <int dim>
- template <int dimm, int spacedim>
+ template <int structdim>
+ template <int dim, int spacedim>
inline
types::global_dof_index
- DoFObjects<dim>::
- get_dof_index (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ DoFObjects<structdim>::
+ get_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index,
const unsigned int fe_index,
const unsigned int local_index,
const unsigned int obj_level) const
{
- Assert ((fe_index != dealii::hp::DoFHandler<dimm,spacedim>::default_fe_index),
+ Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
ExcMessage ("You need to specify a FE index when working "
"with hp DoFHandlers"));
Assert (&dof_handler != 0,
Assert (fe_index < dof_handler.get_fe().size(),
ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
Assert (local_index <
- dof_handler.get_fe()[fe_index].template n_dofs_per_object<dim>(),
+ dof_handler.get_fe()[fe_index].template n_dofs_per_object<structdim>(),
ExcIndexRange(local_index, 0,
dof_handler.get_fe()[fe_index]
- .template n_dofs_per_object<dim>()));
+ .template n_dofs_per_object<structdim>()));
Assert (obj_index < dof_offsets.size(),
ExcIndexRange (obj_index, 0, dof_offsets.size()));
"information for an object on which no such "
"information is available"));
- if (dim == dimm)
- {
- // if we are on a cell, then
- // the only set of indices we
- // store is the one for the
- // cell, which is unique. then
- // fe_index must be
- // active_fe_index
- Assert (fe_index == dof_handler.levels[obj_level]->active_fe_indices[obj_index],
- ExcMessage ("FE index does not match that of the present cell"));
- return dofs[dof_offsets[obj_index]+local_index];
- }
- else
- {
- // we are in higher space
- // dimensions, so there may
- // be multiple finite
- // elements associated with
- // this object. hop along
- // the list of index sets
- // until we find the one
- // with the correct
- // fe_index, and then poke
- // into that part. trigger
- // an exception if we can't
- // find a set for this
- // particular fe_index
- const types::global_dof_index starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- while (true)
- {
- Assert (*pointer != numbers::invalid_dof_index,
- ExcInternalError());
- if (*pointer == fe_index)
- return *(pointer + 1 + local_index);
- else
- pointer += static_cast<types::global_dof_index>(
- dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<dim>() + 1);
- }
- }
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const types::global_dof_index starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ while (true)
+ {
+ Assert (*pointer != numbers::invalid_dof_index,
+ ExcInternalError());
+ if (*pointer == fe_index)
+ return *(pointer + 1 + local_index);
+ else
+ pointer += static_cast<types::global_dof_index>(
+ dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1);
+ }
}
- template <int dim>
- template <int dimm, int spacedim>
+ template <int structdim>
+ template <int dim, int spacedim>
inline
void
- DoFObjects<dim>::
- set_dof_index (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ DoFObjects<structdim>::
+ set_dof_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index,
const unsigned int fe_index,
const unsigned int local_index,
const types::global_dof_index global_index,
const unsigned int obj_level)
{
- Assert ((fe_index != dealii::hp::DoFHandler<dimm,spacedim>::default_fe_index),
+ Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
ExcMessage ("You need to specify a FE index when working "
"with hp DoFHandlers"));
Assert (&dof_handler != 0,
Assert (fe_index < dof_handler.get_fe().size(),
ExcIndexRange (fe_index, 0, dof_handler.get_fe().size()));
Assert (local_index <
- dof_handler.get_fe()[fe_index].template n_dofs_per_object<dim>(),
+ dof_handler.get_fe()[fe_index].template n_dofs_per_object<structdim>(),
ExcIndexRange(local_index, 0,
dof_handler.get_fe()[fe_index]
- .template n_dofs_per_object<dim>()));
+ .template n_dofs_per_object<structdim>()));
Assert (obj_index < dof_offsets.size(),
ExcIndexRange (obj_index, 0, dof_offsets.size()));
"information for an object on which no such "
"information is available"));
- if (dim == dimm)
- {
- // if we are on a cell, then
- // the only set of indices we
- // store is the one for the
- // cell, which is unique. then
- // fe_index must be
- // active_fe_index
- Assert (fe_index == dof_handler.levels[obj_level]->active_fe_indices[obj_index],
- ExcMessage ("FE index does not match that of the present cell"));
- dofs[dof_offsets[obj_index]+local_index] = global_index;
- }
- else
- {
- // we are in higher space
- // dimensions, so there may
- // be multiple finite
- // elements associated with
- // this object. hop along
- // the list of index sets
- // until we find the one
- // with the correct
- // fe_index, and then poke
- // into that part. trigger
- // an exception if we can't
- // find a set for this
- // particular fe_index
- const types::global_dof_index starting_offset = dof_offsets[obj_index];
- types::global_dof_index *pointer = &dofs[starting_offset];
- while (true)
- {
- Assert (*pointer != numbers::invalid_dof_index,
- ExcInternalError());
- if (*pointer == fe_index)
- {
- *(pointer + 1 + local_index) = global_index;
- return;
- }
- else
- pointer += dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<dim>() + 1;
- }
- }
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const types::global_dof_index starting_offset = dof_offsets[obj_index];
+ types::global_dof_index *pointer = &dofs[starting_offset];
+ while (true)
+ {
+ Assert (*pointer != numbers::invalid_dof_index,
+ ExcInternalError());
+ if (*pointer == fe_index)
+ {
+ *(pointer + 1 + local_index) = global_index;
+ return;
+ }
+ else
+ pointer += dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1;
+ }
}
- template <int dim>
- template <int dimm, int spacedim>
+ template <int structdim>
+ template <int dim, int spacedim>
inline
unsigned int
- DoFObjects<dim>::
- n_active_fe_indices (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ DoFObjects<structdim>::
+ n_active_fe_indices (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index) const
{
- Assert (dim <= dimm, ExcInternalError());
Assert (&dof_handler != 0,
ExcMessage ("No DoFHandler is specified for this iterator"));
Assert (&dof_handler.get_fe() != 0,
if (dof_offsets[obj_index] == numbers::invalid_dof_index)
return 0;
- // if we are on a cell, then the
- // only set of indices we store
- // is the one for the cell,
- // which is unique
- if (dim == dimm)
- return 1;
- else
- {
- // otherwise, there may be
- // multiple finite elements
- // associated with this
- // object. hop along the
- // list of index sets until
- // we find the one with the
- // correct fe_index, and
- // then poke into that
- // part. trigger an
- // exception if we can't
- // find a set for this
- // particular fe_index
- const unsigned int starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- unsigned int counter = 0;
- while (true)
- {
- if (*pointer == numbers::invalid_dof_index)
- // end of list reached
- return counter;
- else
- {
- ++counter;
- pointer += dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<dim>() + 1;
- }
- }
- }
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with this
+ // object. hop along the list of index sets until we find the
+ // one with the correct fe_index, and then poke into that
+ // part. trigger an exception if we can't find a set for this
+ // particular fe_index
+ const unsigned int starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ unsigned int counter = 0;
+ while (true)
+ {
+ if (*pointer == numbers::invalid_dof_index)
+ // end of list reached
+ return counter;
+ else
+ {
+ ++counter;
+ pointer += dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1;
+ }
+ }
}
- template <int dim>
- template <int dimm, int spacedim>
+ template <int structdim>
+ template <int dim, int spacedim>
inline
types::global_dof_index
- DoFObjects<dim>::
- nth_active_fe_index (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ DoFObjects<structdim>::
+ nth_active_fe_index (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_level,
const unsigned int obj_index,
const unsigned int n) const
{
- Assert (dim <= dimm, ExcInternalError());
Assert (&dof_handler != 0,
ExcMessage ("No DoFHandler is specified for this iterator"));
Assert (&dof_handler.get_fe() != 0,
"information for an object on which no such "
"information is available"));
- if (dim == dimm)
- {
- // this is a cell, so there
- // is only a single
- // fe_index
- Assert (n == 0, ExcIndexRange (n, 0, 1));
-
- return dof_handler.levels[obj_level]->active_fe_indices[obj_index];
- }
- else
- {
- Assert (n < n_active_fe_indices(dof_handler, obj_index),
- ExcIndexRange (n, 0,
- n_active_fe_indices(dof_handler, obj_index)));
-
- // we are in higher space
- // dimensions, so there may
- // be multiple finite
- // elements associated with
- // this object. hop along
- // the list of index sets
- // until we find the one
- // with the correct
- // fe_index, and then poke
- // into that part. trigger
- // an exception if we can't
- // find a set for this
- // particular fe_index
- const unsigned int starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- unsigned int counter = 0;
- while (true)
- {
- Assert (*pointer != numbers::invalid_dof_index,
- ExcInternalError());
-
- const unsigned int fe_index = *pointer;
-
- Assert (fe_index < dof_handler.get_fe().size(),
- ExcInternalError());
-
- if (counter == n)
- return fe_index;
-
- ++counter;
- pointer += dof_handler.get_fe()[fe_index]
- .template n_dofs_per_object<dim>() + 1;
- }
- }
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ Assert (n < n_active_fe_indices(dof_handler, obj_index),
+ ExcIndexRange (n, 0,
+ n_active_fe_indices(dof_handler, obj_index)));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const unsigned int starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ unsigned int counter = 0;
+ while (true)
+ {
+ Assert (*pointer != numbers::invalid_dof_index,
+ ExcInternalError());
+
+ const unsigned int fe_index = *pointer;
+
+ Assert (fe_index < dof_handler.get_fe().size(),
+ ExcInternalError());
+
+ if (counter == n)
+ return fe_index;
+
+ ++counter;
+ pointer += dof_handler.get_fe()[fe_index]
+ .template n_dofs_per_object<structdim>() + 1;
+ }
}
- template <int dim>
- template <int dimm, int spacedim>
+ template <int structdim>
+ template <int dim, int spacedim>
inline
bool
- DoFObjects<dim>::
- fe_index_is_active (const dealii::hp::DoFHandler<dimm,spacedim> &dof_handler,
+ DoFObjects<structdim>::
+ fe_index_is_active (const dealii::hp::DoFHandler<dim,spacedim> &dof_handler,
const unsigned int obj_index,
const unsigned int fe_index,
const unsigned int obj_level) const
"this DoFHandler"));
Assert (obj_index < dof_offsets.size(),
ExcIndexRange (obj_index, 0, static_cast<unsigned int>(dof_offsets.size())));
- Assert ((fe_index != dealii::hp::DoFHandler<dimm,spacedim>::default_fe_index),
+ Assert ((fe_index != dealii::hp::DoFHandler<dim,spacedim>::default_fe_index),
ExcMessage ("You need to specify a FE index when working "
"with hp DoFHandlers"));
Assert (fe_index < dof_handler.get_fe().size(),
"information for an object on which no such "
"information is available"));
- if (dim == dimm)
- {
- // if we are on a cell,
- // then the only set of
- // indices we store is the
- // one for the cell, which
- // is unique
- Assert (obj_index < dof_handler.levels[obj_level]->active_fe_indices.size(),
- ExcInternalError());
- return (fe_index == dof_handler.levels[obj_level]->active_fe_indices[obj_index]);
- }
- else
- {
- // we are in higher space
- // dimensions, so there may
- // be multiple finite
- // elements associated with
- // this object. hop along
- // the list of index sets
- // until we find the one
- // with the correct
- // fe_index, and then poke
- // into that part. trigger
- // an exception if we can't
- // find a set for this
- // particular fe_index
- const types::global_dof_index starting_offset = dof_offsets[obj_index];
- const types::global_dof_index *pointer = &dofs[starting_offset];
- while (true)
- {
- if (*pointer == numbers::invalid_dof_index)
- // end of list reached
- return false;
- else if (*pointer == fe_index)
- return true;
- else
- pointer += static_cast<types::global_dof_index>(
- dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<dim>()+1);
- }
- }
+ Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
+
+ // there may be multiple finite elements associated with
+ // this object. hop along the list of index sets until we
+ // find the one with the correct fe_index, and then poke
+ // into that part. trigger an exception if we can't find a
+ // set for this particular fe_index
+ const types::global_dof_index starting_offset = dof_offsets[obj_index];
+ const types::global_dof_index *pointer = &dofs[starting_offset];
+ while (true)
+ {
+ if (*pointer == numbers::invalid_dof_index)
+ // end of list reached
+ return false;
+ else if (*pointer == fe_index)
+ return true;
+ else
+ pointer += static_cast<types::global_dof_index>(
+ dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>()+1);
+ }
}
}
*/
reference operator[] (const size_type i);
+ /**
+ * A collective get operation: instead
+ * of getting individual elements of a
+ * vector, this function allows to get
+ * a whole set of elements at once. The
+ * indices of the elements to be read
+ * are stated in the first argument,
+ * the corresponding values are returned in the
+ * second.
+ */
+ template <typename OtherNumber>
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const;
+
+ /**
+ * Just as the above, but with pointers.
+ * Useful in minimizing copying of data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ void extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const;
+
/**
* Copy operator: fill all components of
* the vector with the given scalar
return operator()(i);
}
+
+
+template <typename VectorType>
+template <typename OtherNumber>
+inline
+void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const
+{
+ for (size_type i = 0; i < indices.size(); ++i)
+ values[i] = operator()(indices[i]);
+}
+
+
+
+template <typename VectorType>
+template <typename ForwardIterator, typename OutputIterator>
+inline
+void BlockVectorBase<VectorType>::extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
+{
+ while (indices_begin != indices_end) {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++; values_begin++;
+ }
+}
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
*/
Number &operator [] (const size_type global_index);
+ /**
+ * A collective get operation: instead
+ * of getting individual elements of a
+ * vector, this function allows to get
+ * a whole set of elements at once. The
+ * indices of the elements to be read
+ * are stated in the first argument,
+ * the corresponding values are returned in the
+ * second.
+ */
+ template <typename OtherNumber>
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const;
+
+ /**
+ * Just as the above, but with pointers.
+ * Useful in minimizing copying of data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ void extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const;
+
/**
* Read access to the data field specified by @p local_index. Locally
* owned indices can be accessed with indices
+ template <typename Number>
+ template <typename OtherNumber>
+ inline
+ void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const
+ {
+ for (size_type i = 0; i < indices.size(); ++i)
+ values[i] = operator()(indices[i]);
+ }
+
+
+
+ template <typename Number>
+ template <typename ForwardIterator, typename OutputIterator>
+ inline
+ void Vector<Number>::extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
+ {
+ while (indices_begin != indices_end) {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++; values_begin++;
+ }
+ }
+
+
+
template <typename Number>
inline
Number
void set (const std::vector<size_type> &indices,
const std::vector<PetscScalar> &values);
+ /**
+ * A collective get operation: instead
+ * of getting individual elements of a
+ * vector, this function allows to get
+ * a whole set of elements at once. The
+ * indices of the elements to be read
+ * are stated in the first argument,
+ * the corresponding values are returned in the
+ * second.
+ */
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<PetscScalar> &values) const;
+
+ /**
+ * Just as the above, but with pointers.
+ * Useful in minimizing copying of data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ void extract_subvector_to (const ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const;
+
/**
* A collective add operation: This
* function adds a whole set of values
return comm;
}
+ inline
+ void VectorBase::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<PetscScalar> &values) const
+ {
+ extract_subvector_to(&(indices[0]), &(indices[0]) + indices.size(), &(values[0]));
+ }
+
+ template <typename ForwardIterator, typename OutputIterator>
+ inline
+ void VectorBase::extract_subvector_to (const ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
+ {
+ const PetscInt n_idx = static_cast<PetscInt>(indices_end - indices_begin);
+ if (n_idx == 0)
+ return;
+
+ // if we are dealing
+ // with a parallel vector
+ if (ghosted )
+ {
+
+ int ierr;
+
+ // there is the possibility
+ // that the vector has
+ // ghost elements. in that
+ // case, we first need to
+ // figure out which
+ // elements we own locally,
+ // then get a pointer to
+ // the elements that are
+ // stored here (both the
+ // ones we own as well as
+ // the ghost elements). in
+ // this array, the locally
+ // owned elements come
+ // first followed by the
+ // ghost elements whose
+ // position we can get from
+ // an index set
+ PetscInt begin, end, i;
+ ierr = VecGetOwnershipRange (vector, &begin, &end);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ Vec locally_stored_elements = PETSC_NULL;
+ ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscInt lsize;
+ ierr = VecGetSize(locally_stored_elements, &lsize);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscScalar *ptr;
+ ierr = VecGetArray(locally_stored_elements, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ for (i = 0; i < n_idx; i++) {
+ const unsigned int index = *(indices_begin+i);
+ if ( index>=static_cast<unsigned int>(begin)
+ && index<static_cast<unsigned int>(end) )
+ {
+ //local entry
+ *(values_begin+i) = *(ptr+index-begin);
+ }
+ else
+ {
+ //ghost entry
+ const unsigned int ghostidx
+ = ghost_indices.index_within_set(index);
+
+ Assert(ghostidx+end-begin<(unsigned int)lsize, ExcInternalError());
+ *(values_begin+i) = *(ptr+ghostidx+end-begin);
+ }
+ }
+
+ ierr = VecRestoreArray(locally_stored_elements, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ }
+ // if the vector is local or the
+ // caller, then simply access the
+ // element we are interested in
+ else
+ {
+ int ierr;
+
+ PetscInt begin, end;
+ ierr = VecGetOwnershipRange (vector, &begin, &end);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscScalar *ptr;
+ ierr = VecGetArray(vector, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ for (PetscInt i = 0; i < n_idx; i++) {
+ const unsigned int index = *(indices_begin+i);
+
+ Assert(index>=static_cast<unsigned int>(begin)
+ && index<static_cast<unsigned int>(end), ExcInternalError());
+
+ *(values_begin+i) = *(ptr+index-begin);
+ }
+
+ ierr = VecRestoreArray(vector, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ }
+ }
+
#endif // DOXYGEN
}
indices);
for (TrilinosWrappers::types::int_type i=0; i<n_indices; ++i)
std::cout << indices[i] << " ";
- std::cout << endl << std::endl;
+ std::cout << std::endl << std::endl;
Assert (ierr <= 0,
ExcAccessToNonPresentElement(row, col_index_ptr[0]));
}
TrilinosScalar
operator [] (const size_type index) const;
+ /**
+ * A collective get operation: instead
+ * of getting individual elements of a
+ * vector, this function allows to get
+ * a whole set of elements at once. The
+ * indices of the elements to be read
+ * are stated in the first argument,
+ * the corresponding values are returned in the
+ * second.
+ */
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<TrilinosScalar> &values) const;
+
+ /**
+ * Just as the above, but with pointers.
+ * Useful in minimizing copying of data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ void extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const;
+
/**
* Return the value of the vector
* entry <i>i</i>. Note that this
+ inline
+ void VectorBase::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<TrilinosScalar> &values) const
+ {
+ for (size_type i = 0; i < indices.size(); ++i)
+ values[i] = operator()(indices[i]);
+ }
+
+
+
+ template <typename ForwardIterator, typename OutputIterator>
+ inline
+ void VectorBase::extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
+ {
+ while (indices_begin != indices_end) {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++; values_begin++;
+ }
+ }
+
+
+
inline
VectorBase::iterator
VectorBase::begin()
* Exactly the same as operator().
*/
Number &operator[] (const size_type i);
+
+ /**
+ * A collective get operation: instead
+ * of getting individual elements of a
+ * vector, this function allows to get
+ * a whole set of elements at once. The
+ * indices of the elements to be read
+ * are stated in the first argument,
+ * the corresponding values are returned in the
+ * second.
+ */
+ template <typename OtherNumber>
+ void extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const;
+
+ /**
+ * Just as the above, but with pointers.
+ * Useful in minimizing copying of data around.
+ */
+ template <typename ForwardIterator, typename OutputIterator>
+ void extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const;
//@}
+template <typename Number>
+template <typename OtherNumber>
+inline
+void Vector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
+ std::vector<OtherNumber> &values) const
+{
+ for (size_type i = 0; i < indices.size(); ++i)
+ values[i] = operator()(indices[i]);
+}
+
+
+
+template <typename Number>
+template <typename ForwardIterator, typename OutputIterator>
+inline
+void Vector<Number>::extract_subvector_to (ForwardIterator indices_begin,
+ const ForwardIterator indices_end,
+ OutputIterator values_begin) const
+{
+ while (indices_begin != indices_end) {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++; values_begin++;
+ }
+}
+
+
+
template <typename Number>
inline
Vector<Number> &
template <int dim, int spacedim>
- const std::vector<types::global_dof_index> &
+ const std::vector<unsigned int> &
Triangulation<dim, spacedim>::get_p4est_tree_to_coarse_cell_permutation() const
{
return p4est_tree_to_coarse_cell_permutation;
if (tria)
{
#ifdef DEAL_II_WITH_P4EST
- //this is a distributed Triangulation. We need to traverse the coarse
- //cells in the order p4est does
+ // this is a distributed Triangulation. We need to traverse the coarse
+ // cells in the order p4est does
for (unsigned int c = 0; c < tria->n_cells (0); ++c)
{
- unsigned int coarse_cell_index =
+ const unsigned int coarse_cell_index =
tria->get_p4est_tree_to_coarse_cell_permutation() [c];
const typename DoFHandler<dim>::level_cell_iterator
std::vector<std::vector<unsigned int> >
active_fe_backup(dof_handler.levels.size ());
for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
- std::swap (dof_handler.levels[level]->active_fe_indices, active_fe_backup[level]);
+ std::swap (dof_handler.levels[level]->active_fe_indices,
+ active_fe_backup[level]);
// delete all levels and set them up
// newly, since vectors are
// finite element is used for it
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
- dof_handler.levels[level]->dof_object.dof_offsets
+ dof_handler.levels[level]->dof_offsets
= std::vector<std::vector<types::global_dof_index>::size_type> (
dof_handler.tria->n_raw_lines(level),
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ (std::vector<types::global_dof_index>::size_type)(-1));
types::global_dof_index next_free_dof = 0;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell!=dof_handler.end_active(level); ++cell)
if (!cell->has_children())
{
- dof_handler.levels[level]->dof_object.dof_offsets[cell->index()] = next_free_dof;
+ dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
next_free_dof += cell->get_fe().dofs_per_line;
}
- dof_handler.levels[level]->dof_object.dofs
+ dof_handler.levels[level]->dofs
= std::vector<types::global_dof_index> (next_free_dof,
DoFHandler<dim,spacedim>::invalid_dof_index);
}
if (!cell->has_children())
counter += cell->get_fe().dofs_per_line;
- Assert (dof_handler.levels[level]->dof_object.dofs.size() == counter,
+ Assert (dof_handler.levels[level]->dofs.size() == counter,
ExcInternalError());
Assert (static_cast<unsigned int>
- (std::count (dof_handler.levels[level]->dof_object.dof_offsets.begin(),
- dof_handler.levels[level]->dof_object.dof_offsets.end(),
- DoFHandler<dim,spacedim>::invalid_dof_index))
+ (std::count (dof_handler.levels[level]->dof_offsets.begin(),
+ dof_handler.levels[level]->dof_offsets.end(),
+ (std::vector<types::global_dof_index>::size_type)(-1)))
==
dof_handler.tria->n_raw_lines(level) - dof_handler.tria->n_active_lines(level),
ExcInternalError());
// finite element is used for it
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
- dof_handler.levels[level]->dof_object.dof_offsets
+ dof_handler.levels[level]->dof_offsets
= std::vector<std::vector<types::global_dof_index>::size_type> (
dof_handler.tria->n_raw_quads(level),
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ (std::vector<types::global_dof_index>::size_type)(-1));
types::global_dof_index next_free_dof = 0;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell!=dof_handler.end_active(level); ++cell)
if (!cell->has_children())
{
- dof_handler.levels[level]->dof_object.dof_offsets[cell->index()] = next_free_dof;
+ dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
next_free_dof += cell->get_fe().dofs_per_quad;
}
- dof_handler.levels[level]->dof_object.dofs
+ dof_handler.levels[level]->dofs
= std::vector<types::global_dof_index> (next_free_dof,
DoFHandler<dim,spacedim>::invalid_dof_index);
}
if (!cell->has_children())
counter += cell->get_fe().dofs_per_quad;
- Assert (dof_handler.levels[level]->dof_object.dofs.size() == counter,
+ Assert (dof_handler.levels[level]->dofs.size() == counter,
ExcInternalError());
Assert (static_cast<unsigned int>
- (std::count (dof_handler.levels[level]->dof_object.dof_offsets.begin(),
- dof_handler.levels[level]->dof_object.dof_offsets.end(),
- DoFHandler<dim,spacedim>::invalid_dof_index))
+ (std::count (dof_handler.levels[level]->dof_offsets.begin(),
+ dof_handler.levels[level]->dof_offsets.end(),
+ (std::vector<types::global_dof_index>::size_type)(-1)))
==
dof_handler.tria->n_raw_quads(level) - dof_handler.tria->n_active_quads(level),
ExcInternalError());
dof_handler.faces->lines.dof_offsets
= std::vector<std::vector<types::global_dof_index>::size_type>
(dof_handler.tria->n_raw_lines(),
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ (std::vector<types::global_dof_index>::size_type)(-1));
dof_handler.faces->lines.dofs
= std::vector<types::global_dof_index> (n_line_slots,
DoFHandler<dim,spacedim>::invalid_dof_index);
// finite element is used for it
for (unsigned int level=0; level<dof_handler.tria->n_levels(); ++level)
{
- dof_handler.levels[level]->dof_object.dof_offsets
+ dof_handler.levels[level]->dof_offsets
= std::vector<std::vector<types::global_dof_index>::size_type>
(dof_handler.tria->n_raw_hexs(level),
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ (std::vector<types::global_dof_index>::size_type)(-1));
types::global_dof_index next_free_dof = 0;
for (typename DoFHandler<dim,spacedim>::active_cell_iterator
cell!=dof_handler.end_active(level); ++cell)
if (!cell->has_children())
{
- dof_handler.levels[level]->dof_object.dof_offsets[cell->index()] = next_free_dof;
+ dof_handler.levels[level]->dof_offsets[cell->index()] = next_free_dof;
next_free_dof += cell->get_fe().dofs_per_hex;
}
- dof_handler.levels[level]->dof_object.dofs
+ dof_handler.levels[level]->dofs
= std::vector<types::global_dof_index> (next_free_dof,
DoFHandler<dim,spacedim>::invalid_dof_index);
}
if (!cell->has_children())
counter += cell->get_fe().dofs_per_hex;
- Assert (dof_handler.levels[level]->dof_object.dofs.size() == counter,
+ Assert (dof_handler.levels[level]->dofs.size() == counter,
ExcInternalError());
Assert (static_cast<unsigned int>
- (std::count (dof_handler.levels[level]->dof_object.dof_offsets.begin(),
- dof_handler.levels[level]->dof_object.dof_offsets.end(),
- DoFHandler<dim,spacedim>::invalid_dof_index))
+ (std::count (dof_handler.levels[level]->dof_offsets.begin(),
+ dof_handler.levels[level]->dof_offsets.end(),
+ (std::vector<types::global_dof_index>::size_type)(-1)))
==
dof_handler.tria->n_raw_hexs(level) - dof_handler.tria->n_active_hexs(level),
ExcInternalError());
dof_handler.faces->quads.dof_offsets
= std::vector<std::vector<types::global_dof_index>::size_type>
(dof_handler.tria->n_raw_quads(),
- DoFHandler<dim,spacedim>::invalid_dof_index);
+ (std::vector<types::global_dof_index>::size_type)(-1));
dof_handler.faces->quads.dofs
= std::vector<types::global_dof_index> (n_quad_slots,
DoFHandler<dim,spacedim>::invalid_dof_index);
std::size_t
DoFLevel<dim>::memory_consumption () const
{
- return MemoryConsumption::memory_consumption (active_fe_indices) +
- MemoryConsumption::memory_consumption (dof_object);
+ return (MemoryConsumption::memory_consumption (active_fe_indices) +
+ MemoryConsumption::memory_consumption (dofs) +
+ MemoryConsumption::memory_consumption (dof_offsets));
}
for (DH : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS)
{
- template class DataOut_DoFData<DH<deal_II_dimension>,deal_II_dimension>;
- template class DataOut_DoFData<DH<deal_II_dimension>,deal_II_dimension+1>;
-
- #if deal_II_dimension < 3
- template class DataOut_DoFData<DH<deal_II_dimension>,deal_II_dimension,deal_II_dimension+1>;
- template class DataOut_DoFData<DH<deal_II_dimension,deal_II_dimension+1>,deal_II_dimension,deal_II_dimension+1>;
- #endif
-
-
template class DataOut<deal_II_dimension, DH<deal_II_dimension> >;
#if deal_II_dimension < 3
template class DataOut<deal_II_dimension, DH<deal_II_dimension,deal_II_dimension+1> >;
projects / dealii-svn.git / commitdiff