#include <deal.II/lac/sparse_direct.h>
// This includes the library for the incomplete LU factorization that will be
-// used as a preconditioner in 3D:
+// used as a preconditioner in 3d:
#include <deal.II/lac/sparse_ilu.h>
// This is C++:
template <int dim>
struct InnerPreconditioner;
- // In 2D, we are going to use a sparse direct solver as preconditioner:
+ // In 2d, we are going to use a sparse direct solver as preconditioner:
template <>
struct InnerPreconditioner<2>
{
using type = SparseDirectUMFPACK;
};
- // And the ILU preconditioning in 3D, called by SparseILU:
+ // And the ILU preconditioning in 3d, called by SparseILU:
template <>
struct InnerPreconditioner<3>
{
// pattern objects.
//
// We then proceed with distributing degrees of freedom and renumbering
- // them: In order to make the ILU preconditioner (in 3D) work efficiently,
+ // them: In order to make the ILU preconditioner (in 3d) work efficiently,
// it is important to enumerate the degrees of freedom in such a way that it
// reduces the bandwidth of the matrix, or maybe more importantly: in such a
// way that the ILU is as close as possible to a real LU decomposition. On
// the same way as in step-20, i.e. directly build an object of type
// SparsityPattern through DoFTools::make_sparsity_pattern. However, there
// is a major reason not to do so:
- // In 3D, the function DoFTools::max_couplings_between_dofs yields a
+ // In 3d, the function DoFTools::max_couplings_between_dofs yields a
// conservative but rather large number for the coupling between the
// individual dofs, so that the memory initially provided for the creation
// of the sparsity pattern of the matrix is far too much -- so much actually
// that the initial sparsity pattern won't even fit into the physical memory
- // of most systems already for moderately-sized 3D problems, see also the
+ // of most systems already for moderately-sized 3d problems, see also the
// discussion in step-18. Instead, we first build temporary objects that use
// a different data structure that doesn't require allocating more memory
// than necessary but isn't suitable for use as a basis of SparseMatrix or
// We then apply an initial refinement before solving for the first
- // time. In 3D, there are going to be more degrees of freedom, so we
+ // time. In 3d, there are going to be more degrees of freedom, so we
// refine less there:
triangulation.refine_global(4 - dim);
// @sect3{Initial conditions}
// In the following two classes, we first implement the exact solution for
- // 1D, 2D, and 3D mentioned in the introduction to this program. This
+ // 1D, 2d, and 3d mentioned in the introduction to this program. This
// space-time solution may be of independent interest if one wanted to test
// the accuracy of the program by comparing the numerical against the
// analytic solution (note however that the program uses a finite domain,
// Next, two points are defined for position and focal point of the
// transducer lens, which is the center of the circle whose segment will
// form the transducer part of the boundary. Notice that this is the only
- // point in the program where things are slightly different in 2D and 3D.
- // Even though this tutorial only deals with the 2D case, the necessary
- // additions to make this program functional in 3D are so minimal that we
+ // point in the program where things are slightly different in 2d and 3d.
+ // Even though this tutorial only deals with the 2d case, the necessary
+ // additions to make this program functional in 3d are so minimal that we
// opt for including them:
const Point<dim> transducer =
(dim == 2) ? Point<dim>(0.5, 0.0) : Point<dim>(0.5, 0.5, 0.0);
face->set_manifold_id(1);
}
// For the circle part of the transducer lens, a SphericalManifold object
- // is used (which, of course, in 2D just represents a circle), with center
+ // is used (which, of course, in 2d just represents a circle), with center
// computed as above.
triangulation.set_manifold(1, SphericalManifold<dim>(focal_point));
// with an iterative solver and a preconditioner that do a good job on this
// matrix. We chose instead to go a different way and solve the linear
// system with the sparse LU decomposition provided by UMFPACK. This is
- // often a good first choice for 2D problems and works reasonably well even
+ // often a good first choice for 2d problems and works reasonably well even
// for a large number of DoFs. The deal.II interface to UMFPACK is given by
// the SparseDirectUMFPACK class, which is very easy to use and allows us to
// solve our linear system with just 3 lines of code.
// Ok, let's start: we need a quadrature formula for the evaluation of the
// integrals on each cell. Let's take a Gauss formula with two quadrature
// points in each direction, i.e. a total of four points since we are in
- // 2D. This quadrature formula integrates polynomials of degrees up to three
- // exactly (in 1D). It is easy to check that this is sufficient for the
+ // 2d. This quadrature formula integrates polynomials of degrees up to three
+ // exactly (in 1d). It is easy to check that this is sufficient for the
// present problem:
QGauss<2> quadrature_formula(fe.degree + 1);
// And we initialize the object which we have briefly talked about above. It
// For use further down below, we define a shortcut for a value that will
// be used very frequently. Namely, an abbreviation for the number of degrees
- // of freedom on each cell (since we are in 2D and degrees of freedom are
+ // of freedom on each cell (since we are in 2d and degrees of freedom are
// associated with vertices only, this number is four, but we rather want to
// write the definition of this variable in a way that does not preclude us
// from later choosing a different finite element that has a different
{
using namespace Step30;
- // If you want to run the program in 3D, simply change the following
+ // If you want to run the program in 3d, simply change the following
// line to <code>const unsigned int dim = 3;</code>.
const unsigned int dim = 2;
// Now let's turn to the temperature part: First, we compute the time step
- // size. We found that we need smaller time steps for 3D than for 2D for
+ // size. We found that we need smaller time steps for 3d than for 2d for
// the shell geometry. This is because the cells are more distorted in
// that case (it is the smallest edge length that determines the CFL
// number). Instead of computing the time step from maximum velocity and
// @sect4{LaplaceProblem::run}
// The function that runs the program is very similar to the one in
- // step-16. We do few refinement steps in 3D compared to 2D, but that's
+ // step-16. We do few refinement steps in 3d compared to 2d, but that's
// it.
//
// Before we run the program, we output some information about the detected
//
// But back to the concrete case here:
// For this tutorial, we choose as right hand side the function
-// $4(x^4+y^4)$ in 2D, or $4(x^4+y^4+z^4)$ in 3D. We could write this
+// $4(x^4+y^4)$ in 2d, or $4(x^4+y^4+z^4)$ in 3d. We could write this
// distinction using an if-statement on the space dimension, but here is a
-// simple way that also allows us to use the same function in 1D (or in 4D, if
+// simple way that also allows us to use the same function in 1d (or in 4D, if
// you should desire to do so), by using a short loop. Fortunately, the
// compiler knows the size of the loop at compile time (remember that at the
// time when you define the template, the compiler doesn't know the value of
}
-// As boundary values, we choose $x^2+y^2$ in 2D, and $x^2+y^2+z^2$ in 3D. This
+// As boundary values, we choose $x^2+y^2$ in 2d, and $x^2+y^2+z^2$ in 3d. This
// happens to be equal to the square of the vector from the origin to the
// point at which we would like to evaluate the function, irrespective of the
// dimension. So that is what we return:
// @sect4{Step4::make_grid}
// Grid creation is something inherently dimension dependent. However, as long
-// as the domains are sufficiently similar in 2D or 3D, the library can
+// as the domains are sufficiently similar in 2d or 3d, the library can
// abstract for you. In our case, we would like to again solve on the square
-// $[-1,1]\times [-1,1]$ in 2D, or on the cube $[-1,1] \times [-1,1] \times
-// [-1,1]$ in 3D; both can be termed GridGenerator::hyper_cube(), so we may
+// $[-1,1]\times [-1,1]$ in 2d, or on the cube $[-1,1] \times [-1,1] \times
+// [-1,1]$ in 3d; both can be termed GridGenerator::hyper_cube(), so we may
// use the same function in whatever dimension we are. Of course, the
// functions that create a hypercube in two and three dimensions are very much
// different, but that is something you need not care about. Let the library
// Next, we again have to loop over all cells and assemble local
// contributions. Note, that a cell is a quadrilateral in two space
- // dimensions, but a hexahedron in 3D. In fact, the
+ // dimensions, but a hexahedron in 3d. In fact, the
// <code>active_cell_iterator</code> data type is something different,
// depending on the dimension we are in, but to the outside world they look
// alike and you will probably never see a difference. In any case, the real
// @sect4{ObstacleProblem::make_grid}
// We solve our obstacle problem on the square $[-1,1]\times [-1,1]$ in
- // 2D. This function therefore just sets up one of the simplest possible
+ // 2d. This function therefore just sets up one of the simplest possible
// meshes.
template <int dim>
void ObstacleProblem<dim>::make_grid()
// in displacement is non-constant between each time step.
constraints.clear();
- // The boundary conditions for the indentation problem in 3D are as
+ // The boundary conditions for the indentation problem in 3d are as
// follows: On the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry
// condition to allow only planar movement while the +x and +z faces
// (IDs 1,5) are traction free. In this contrived problem, part of the
// We define a time-dependent function that is used as initial
// value. Different solutions can be obtained by varying the starting
// time. This function, taken from step-25, would represent an analytic
- // solution in 1D for all times, but is merely used for setting some
+ // solution in 1d for all times, but is merely used for setting some
// starting solution of interest here. More elaborate choices that could
// test the convergence of this program are given in step-25.
template <int dim>
// two-dimensional triangulation, while this function is a template for
// arbitrary dimension. Since this is only a demonstration program, we will
// not use different input files for the different dimensions, but rather
- // quickly kill the whole program if we are not in 2D. Of course, since the
+ // quickly kill the whole program if we are not in 2d. Of course, since the
// main function below assumes that we are working in two dimensions we
// could skip this check, in this version of the program, without any ill
// effects.
// Next comes the implementation of the convection velocity. As described in
- // the introduction, we choose a velocity field that is $(y, -x)$ in 2D and
- // $(y, -x, 1)$ in 3D. This gives a divergence-free velocity field.
+ // the introduction, we choose a velocity field that is $(y, -x)$ in 2d and
+ // $(y, -x, 1)$ in 3d. This gives a divergence-free velocity field.
template <int dim>
class ConvectionVelocity : public TensorFunction<1, dim>
{
// polynomial degree $k$ in $d$ space dimensions, out of the $(k+1)^d$
// degrees of freedom of a cell. A similar reduction is also possible for
// the interior penalty method that evaluates values and first derivatives
- // on the faces. When using a Hermite-like basis in 1D, only up to two basis
+ // on the faces. When using a Hermite-like basis in 1d, only up to two basis
// functions contribute to the value and derivative. The class FE_DGQHermite
// implements a tensor product of this concept, as discussed in the
// introduction. Thus, only $2(k+1)^{d-1}$ degrees of freedom must be
// Next we turn to the preconditioner initialization. As explained in the
// introduction, we want to construct an (approximate) inverse of the cell
- // matrices from a product of 1D mass and Laplace matrices. Our first task
- // is to compute the 1D matrices, which we do by first creating a 1D finite
+ // matrices from a product of 1d mass and Laplace matrices. Our first task
+ // is to compute the 1d matrices, which we do by first creating a 1d finite
// element. Instead of anticipating FE_DGQHermite<1> here, we get the finite
// element's name from DoFHandler, replace the @p dim argument (2 or 3) by 1
- // to create a 1D name, and construct the 1D element by using FETools.
+ // to create a 1d name, and construct the 1d element by using FETools.
template <int dim, int fe_degree, typename number>
void PreconditionBlockJacobi<dim, fe_degree, number>::initialize(
name.replace(name.find('<') + 1, 1, "1");
std::unique_ptr<FiniteElement<1>> fe_1d = FETools::get_fe_by_name<1>(name);
- // As for computing the 1D matrices on the unit element, we simply write
+ // As for computing the 1d matrices on the unit element, we simply write
// down what a typical assembly procedure over rows and columns of the
// matrix as well as the quadrature points would do. We select the same
// Laplace matrices once and for all using the coefficients 0.5 for
// The left and right boundary terms assembled by the next two
// statements appear to have somewhat arbitrary signs, but those are
// correct as can be verified by looking at step-39 and inserting
- // the value -1 and 1 for the normal vector in the 1D case.
+ // the value -1 and 1 for the normal vector in the 1d case.
sum_laplace +=
(1. * fe_1d->shape_value(i, Point<1>()) *
fe_1d->shape_value(j, Point<1>()) * op.get_penalty_factor() +
// we check that it is diagonal and then extract the determinant of the
// original Jacobian, i.e., the inverse of the determinant of the inverse
// Jacobian, and set the weight as $\text{det}(J) / h_d^2$ according to
- // the 1D Laplacian times $d-1$ copies of the mass matrix.
+ // the 1d Laplacian times $d-1$ copies of the mass matrix.
cell_matrices.clear();
FEEvaluation<dim, fe_degree, fe_degree + 1, 1, number> phi(*data);
unsigned int old_mapping_data_index = numbers::invalid_unsigned_int;
// periodic boundary conditions in the $x$-direction, a Dirichlet condition
// on the front face in $y$ direction (i.e., the face with index number 2,
// with boundary id equal to 0), and Neumann conditions on the back face as
- // well as the two faces in $z$ direction for the 3D case (with boundary id
+ // well as the two faces in $z$ direction for the 3d case (with boundary id
// equal to 1). The extent of the domain is a bit different in the $x$
// direction (where we want to achieve a periodic solution given the
// definition of `Solution`) as compared to the $y$ and $z$ directions.
// for degree 1 and degree 3 finite elements. If the user wants to change to
// another degree, they may need to adjust these numbers. For block smoothers,
// this parameter has a more straightforward interpretation, namely that for
- // additive methods in 2D, a DoF can have a repeated contribution from up to 4
+ // additive methods in 2d, a DoF can have a repeated contribution from up to 4
// cells, therefore we must relax these methods by 0.25 to compensate. This is
// not an issue for multiplicative methods as each cell's inverse application
// carries new information to all its DoFs.
// ball should be compatible with the outer grid in the sense that their
// vertices coincide so as to allow the two grid to be merged. The grid coming
// out of GridGenerator::hyper_shell fulfills the requirements on the inner
- // side in case it is created with $2d$ coarse cells (6 coarse cells in 3D
+ // side in case it is created with $2d$ coarse cells (6 coarse cells in 3d
// which we are going to use) – this is the same number of cells as
// there are boundary faces for the ball. For the outer surface, we use the
// fact that the 6 faces on the surface of the shell without a manifold
//
// The reason for choosing this somewhat unusual scheme is due to the heavy
// work involved in computing the cell matrix for a relatively high
- // polynomial degree in 3D. As we want to highlight the cost of the mapping
+ // polynomial degree in 3d. As we want to highlight the cost of the mapping
// in this tutorial program, we better do the assembly in an optimized way
// in order to not chase bottlenecks that have been solved by the community
// already. Matrix-matrix multiplication is one of the best optimized
// For the parallel computation we define a
// parallel::distributed::Triangulation. As the computational domain is a
- // circle in 2D and a ball in 3D, we assign in addition to the
+ // circle in 2d and a ball in 3d, we assign in addition to the
// SphericalManifold for boundary cells a TransfiniteInterpolationManifold
// object for the mapping of the inner cells, which takes care of the inner
// cells. In this example we use an isoparametric finite element approach
// framework and disable shared-memory parallelization by limiting the number of
// threads to one. Finally to run the solver for the <i>Gelfand problem</i> we
// create an object of the <code>GelfandProblem</code> class and call the run
-// function. Exemplarily we solve the problem once in 2D and once in 3D each
+// function. Exemplarily we solve the problem once in 2d and once in 3d each
// with fourth-order Lagrangian finite elements.
int main(int argc, char *argv[])
{
// want to use for the nonlinear terms in the Euler equations. Furthermore,
// we specify the time interval for the time-dependent problem, and
// implement two different test cases. The first one is an analytical
- // solution in 2D, whereas the second is a channel flow around a cylinder as
+ // solution in 2d, whereas the second is a channel flow around a cylinder as
// described in the introduction. Depending on the test case, we also change
// the final time up to which we run the simulation, and a variable
// `output_tick` that specifies in which intervals we want to write output
// step-33. The interface of the DataPostprocessor class is intuitive,
// requiring us to provide information about what needs to be evaluated
// (typically only the values of the solution, except for the Schlieren plot
- // that we only enable in 2D where it makes sense), and the names of what
+ // that we only enable in 2d where it makes sense), and the names of what
// gets evaluated. Note that it would also be possible to extract most
// information by calculator tools within visualization programs such as
// ParaView, but it is so much more convenient to do it already when writing
// choose a constant inflow profile, whereas we set a subsonic outflow at
// the right. For the boundary around the cylinder (boundary id equal to 2)
// as well as the channel walls (boundary id equal to 3) we use the wall
- // boundary type, which is no-normal flow. Furthermore, for the 3D cylinder
+ // boundary type, which is no-normal flow. Furthermore, for the 3d cylinder
// we also add a gravity force in vertical direction. Having the base mesh
// in place (including the manifolds set by
// GridGenerator::channel_with_cylinder()), we can then perform the
//
// For the purpose of this example step
// we simply implement a homogeneous uniform flow field for which the
- // direction and a 1D primitive state (density, velocity, pressure) are
+ // direction and a 1d primitive state (density, velocity, pressure) are
// read from the parameter file.
//
// It would be desirable to initialize the class in a single shot:
// We start again by defining a couple of helper functions:
//
// The first function takes a state <code>U</code> and a unit vector
- // <code>n_ij</code> and computes the <i>projected</i> 1D state in
+ // <code>n_ij</code> and computes the <i>projected</i> 1d state in
// direction of the unit vector.
namespace
{
const auto perpendicular_m = m - projected_U[1] * n_ij;
projected_U[2] = U[1 + dim] - 0.5 * perpendicular_m.norm_square() / U[0];
- // We return the 1D state in <i>primitive</i> variables instead of
+ // We return the 1d state in <i>primitive</i> variables instead of
// conserved quantities. The return array consists of density $\rho$,
// velocity $u$, pressure $p$ and local speed of sound $a$:
// <code>1</code>, which is why you haven't seen this parameter in
// previous examples). This parameter denotes into how many sub-cells per
// space direction each cell shall be subdivided for output. For example,
- // if you give <code>2</code>, this leads to 4 cells in 2D and 8 cells in
- // 3D. For quadratic elements, two sub-cells per space direction is
+ // if you give <code>2</code>, this leads to 4 cells in 2d and 8 cells in
+ // 3d. For quadratic elements, two sub-cells per space direction is
// obviously the right choice, so this is what we choose. In general, for
// elements of polynomial order <code>q</code>, we use <code>q</code>
// subdivisions, and the order of the elements is determined in the same
"extension. Bailing out."));
// Now we check how many faces are contained in the `Shape`. OpenCASCADE
- // is intrinsically 3D, so if this number is zero, we interpret this as
+ // is intrinsically 3d, so if this number is zero, we interpret this as
// a line manifold, otherwise as a
// OpenCASCADE::NormalToMeshProjectionManifold in `spacedim` = 3, or
// OpenCASCADE::NURBSPatchManifold in `spacedim` = 2.
// We will use the FEEvaluation class to evaluate the solution vector
// at the quadrature points and to perform the integration. In contrast to
// other tutorials, the template arguments `degree` is set to $-1$ and
- // `number of quadrature in 1D` to $0$. In this case, FEEvaluation selects
+ // `number of quadrature in 1d` to $0$. In this case, FEEvaluation selects
// dynamically the correct polynomial degree and number of quadrature
// points. Here, we introduce an alias to FEEvaluation with the correct
// template parameters so that we do not have to worry about them later on.
// @sect4{LaplaceProblem::initialize_grid}
// For a L-shaped domain, we could use the function GridGenerator::hyper_L()
- // as demonstrated in step-50. However in the 2D case, that particular
+ // as demonstrated in step-50. However in the 2d case, that particular
// function removes the first quadrant, while we need the fourth quadrant
// removed in our scenario. Thus, we will use a different function
// GridGenerator::subdivided_hyper_L() which gives us more options to create
// the mesh. Furthermore, we formulate that function in a way that it also
- // generates a 3D mesh: the 2D L-shaped domain will basically elongated by 1
+ // generates a 3d mesh: the 2d L-shaped domain will basically elongated by 1
// in the positive z-direction.
//
// We first pretend to build a GridGenerator::subdivided_hyper_rectangle().
phi_m.reinit(cell, face);
// Interpolate the values from the cell quadrature points to the
- // quadrature points of the current face via a simple 1D
+ // quadrature points of the current face via a simple 1d
// interpolation:
internal::FEFaceNormalEvaluationImpl<dim,
n_points_1d - 1,
}
// Evaluate local integrals related to cell by quadrature and
- // add into cell contribution via a simple 1D interpolation:
+ // add into cell contribution via a simple 1d interpolation:
internal::FEFaceNormalEvaluationImpl<dim,
n_points_1d - 1,
VectorizedArrayType>::
// The reason for refining is a bit accidental: we use the QGauss
// quadrature formula with two points in each direction for integration of the
// right hand side; that means that there are four quadrature points on each
- // cell (in 2D). If we only refine the initial grid once globally, then there
+ // cell (in 2d). If we only refine the initial grid once globally, then there
// will be only four quadrature points in each direction on the
// domain. However, the right hand side function was chosen to be rather
// localized and in that case, by pure chance, it happens that all quadrature
// This is the <code>main</code> function. We define here the number of mesh
// refinements, the polynomial degree for the two finite element spaces
// (for the solution and the two liftings) and the two penalty coefficients.
-// We can also change the dimension to run the code in 3D.
+// We can also change the dimension to run the code in 3d.
int main()
{
try
// To make the assembly easier, we use the class NonMatching::FEValues,
// which does the above steps 1 and 2 for us. The algorithm @cite saye_2015
// that is used to generate the quadrature rules on the intersected cells
- // uses a 1-dimensional quadrature rule as base. Thus, we pass a 1D
+ // uses a 1-dimensional quadrature rule as base. Thus, we pass a 1d
// Gauss--Legendre quadrature to the constructor of NonMatching::FEValues.
- // On the non-intersected cells, a tensor product of this 1D-quadrature will
+ // On the non-intersected cells, a tensor product of this 1d-quadrature will
// be used.
//
// As stated in the introduction, each cell has 3 different regions: inside,
// understand data that is associated with nodes: they cannot plot
// fifth-degree basis functions, which results in a very inaccurate picture
// of the solution we computed. To get around this we save multiple
- // <em>patches</em> per cell: in 2D we save 64 bilinear `cells' to the VTU
- // file for each cell, and in 3D we save 512. The end result is that the
+ // <em>patches</em> per cell: in 2d we save 64 bilinear `cells' to the VTU
+ // file for each cell, and in 3d we save 512. The end result is that the
// visualization program will use a piecewise linear interpolation of the
// cubic basis functions: this captures the solution detail and, with most
// screen resolutions, looks smooth. We save the grid in a separate step
// Then check whether the neighbor is active. If it is, then it
// is on the same level or one level coarser (if we are not in
- // 1D), and we are interested in it in any case.
+ // 1d), and we are interested in it in any case.
if (neighbor->is_active())
scratch_data.active_neighbors.push_back(neighbor);
else
projects / dealii.git / commitdiff