# Finite Element Programming

## Introduction

NDSolve provides a high-level, one-step interface for solving partial differential equations with the finite element method. However, you may want to control the steps of the solution process with more detail. The NDSolve`FEM` package provides a lower-level interface that gives extensive control for each part of the solution process.

In[1]:= |

The low-level functions in the NDSolve`FEM` package may be used for a variety of purposes:

- to enable development of specific, finite-element-based solvers, not only to solve PDEs but also other areas of numerics

- to use NDSolve as an equation preprocessor

## Finite Element Data within NDSolve

The NDSolve framework has lower-level functions NDSolve`ProcessEquations, NDSolve`Iterate, and NDSolve`ProcessSolutions and a data object NDSolve`StateData that provide some access to the solution process for all types of differential equations. The general use of these functions is described in "Components and Data Structures".

Consistent with the rest of NDSolve, the finite element method has its data accessible through NDSolve`StateData. The following section explores how to access the finite element data from the NDSolve`StateData object. The finite element data is stored in an NDSolve`FEM`FiniteElementData object and accessed through the "FiniteElementData" property.

When there is no time dependence in the problem, as in this case, the display form of the NDSolve`StateData object will indicate this by displaying "SteadyState". For transient equations, the time interval over which it is to be integrated will be shown.

The finite element data can be extracted from the state object.

All data objects generated during the solution process are documented, for example, FiniteElementData, as are their constructor functions. This tutorial will explain the usage of these functions.

The display form for a FiniteElementData object shows the degrees of freedom of the discretized partial differential equation (PDE). The degrees of freedom are typically the number of nodes in the mesh times the number of dependent variables. The degrees of freedom indicate the number of rows and columns of the matrices that are to be assembled.

For a steady state problem, invoking NDSolve`Iterate finds the solution of a system of equations with LinearSolve.

In[3]:= |

The solution is then stored in the finite element data object.

The solution is represented as an ×1 matrix, where is the degrees of freedom.

## Passing Finite Element Options to NDSolve

In the subsequent sections, the finite element functions internal to NDSolve will be discussed. Some of these functions have options. It is always possible to pass these options to NDSolve without making use of internal finite element functions. These are given in the Method option to NDSolve.

NDSolve[…,Method{"PDEDiscretization"{"FiniteElement", femopts}}] | use the finite element method with options femopts for a steady-state problem |

NDSolve[…,Method{"PDEDiscretization"{"MethodOfLines", "SpatialDiscretization"->{"FiniteElement", femopts}}}] | use the finite element method with options femopts for the spatial discretization of a transient problem |

Specifying finite element method options.

The same is possible for transient PDEs.

## A Workflow Overview

The following sections show how NDSolve solves a finite element model step by step. NDSolve essentially uses the same functions as seen in these sections. Furthermore, the content of the NDSolve`FEM`FiniteElementData object is created through the functions introduced below.

The solution is found in three stages:

During the initialization stage, the PDE and boundary conditions are analyzed and classified into different components, and these results are stored in the PDECoefficientData and BoundaryConditionData objects, respectively. The finite element method data is set up and stored in the FEMMethodData object.

InitializePDECoefficients | initialize partial differential equation coefficients |

InitializeBoundaryConditions | initialize boundary conditions |

InitializePDEMethodData | initialize partial differential equation method data |

In the second stage, the PDE and boundary conditions are discretized and stored in DiscretizedPDEData and DiscretizedBoundaryConditionData. "Discretized" essentially means that the continuous PDE and boundary conditions are approximated by discrete versions represented by so-called system matrices.

DiscretizePDE | discretize initialized partial differential equations |

DiscretizeBoundaryConditions | discretize initialized boundary conditions |

The last step is then to merge the discretized PDE with the discretized boundary conditions (DeployBoundaryConditions) and to use LinearSolve to solve the equations.

DeployBoundaryConditions | deploy discretized boundary conditions into discretized partial differential equations |

LinearSolve | solve linear systems of equations |

Each of these stages will be addressed. Before that, however, a model equation is set up.

#### The Partial Differential Equation Problem Setup

To make a PDE susceptible to be solved by a numerical method such as the finite element method, three components are needed:

The following setup will be used as a model PDE in the subsequent explanations.

To solve the PDE over the region with the given boundary conditions, NDSolve can be used directly.

This is a typical workflow that uses the one-step capability of NDSolve.

To compare with using the NDSolve`FEM` functions to solve step by step, start by extracting and inspecting the data used by NDSolve.

The PDECoefficientData has been created by a call to InitializePDECoefficients. The BoundaryConditionData has been created by a call to InitializeBoundaryConditions, and the FEMMethodData has been created by a call to InitializePDEMethodData. These data objects hold data for subsequent computations. The creation of these data objects is discussed in the next section.

## Stationary PDEs

### Initialization Stage

The initialization stage creates all data objects that are created during a run of NDSolve`ProcessEquations. This section will illustrate that process. It is possible to skip this section and continue with the discretization stage and make use of the initialized data structures NDSolve`ProcessEquations creates. With this it is possible to use NDSolve`ProcessEquations as an equation preprocessor, for example, for a new numerical discretization method.

##### Variable and Solution Data

For NDSolve to find a solution, the specification of the dependent (u) and independent ({x,y}) variables and the region (Ω) needs to be given. The same is true for the finite element functions. The variable and solution data are stored in lists that make a specific data structure corresponding to different components. These data lists may be generated using NDSolve`VariableData and NDSolve`SolutionData. Elements of the data lists can be filled in or modified using NDSolve`SetSolutionDataComponent.

Accompanying the variable data is a solution data structure. The solution data structure is essentially the numerical incarnation of the variable data. For example, the numerical version of the "Space" component of the solution data is the NumericalRegion that will also contain the ElementMesh. "Element Mesh Generation" gives more information about the creation of the mesh from the region Ω.

##### InitializePDECoefficients

Consider the following partial differential equation in one dependent variable :

That PDE is defined in . The coefficients , , , and are scalars. , , and are vectors of length and is an × matrix. The input coefficients to InitializePDECoefficients are found by correspondence of the above PDE with the PDE to be modeled.

All coefficients that are not explicitly specified are assumed to be zero.

The initialized coefficients are stored in the PDECoefficientData data structure. The display form shows the system size, i.e. the number of dependent variables (1), and the operator size, i.e the space dimension (2) of the initialized coefficients.

During the initialization, the coefficients are classified into different categories of coefficients, such as stationary coefficients and transient coefficients. This will be discussed later in more detail.

##### InitializeBoundaryConditions

Similar to the partial differential equation initialization, the boundary conditions need to be initialized.

Boundary condition coefficients are classified into categories similar to the PDE coefficients.

##### InitializePDEMethodData

To initialize the finite element method data, which is needed in subsequent discretization steps, InitializePDEMethodData is used. Currently, the only discretization method available in this framework is the finite element method. Thus, by default, InitializePDEMethodData generates a FEMMethodData object. During the initialization, a method-specific setup is performed; for example, among other things, the interpolation and integration order of the method are set up and an ElementMesh is generated from the NumericalRegion.

The initialized PDE method data is stored in the FEMMethodData object. The display form of FEMMethodData shows degrees of freedom, the interpolation order for each dependent variable, and integration order.

All mesh generation options as discussed in the "Element Mesh Generation" tutorial can be specified.

During the method initialization, an ElementMesh object will be generated and stored in the NumericalRegion object.

The degrees of freedom of a system of equations result from a combination of the nodes in the element mesh, the number of dependent variables, and the interpolation order of the dependent variables if multiple dependent variables are given. When one dependent variable is given, the degrees of freedom correspond to the nodes in the ElementMesh object.

The "IntegrationOrder" is the order of accuracy used to integrate the finite element operators.

For a first-order element mesh, the default integration order is to use second-order accurate integration points, while for a second-order element mesh, the default is to use fourth-order accurate integration points. The maximum supported integration order is fifth order.

Reducing the integration order will result in a faster solution, but potentially a less-accurate solution. Increasing the integration order will result in a longer time to solve the finite element model.

In[35]:= |

##### Extracting the Initialized Finite Element Data from the NDSolve`StateData

As an alternative to the manual initialization of the finite element method data, PDE coefficients, and boundary conditions, NDSolve`ProcessEquations can be utilized for this purpose.

In some cases, details like arithmetic order or generated unique variable names may be different, but otherwise the expressions are the same.

### Discretization Stage

During the discretization stage, the continuous PDE and boundary conditions are approximated with discrete variants. This process is at the heart of the partial differential equation solving process.

##### DiscretizePDE

The initialized PDE coefficients are transformed into a discrete form.

The display form of the DiscretizedPDEData shows the first dimension of the system matrices.

During the matrix assembly process, access to the finite elements is possible. Note that the form of the finite elements is such that the matrix assembly process is efficient.

##### DiscretizeBoundaryCondition

DiscretizeBoundaryConditions computes a discretized version of the boundary conditions given. In a later step, the discretized boundary conditions are then deployed into the actual system matrices.

The DirichletCondition and NeumannValue boundary specifications are represented differently. The contributions from NeumannValue are to the system load matrix and the system stiffness matrix.

The "StiffnessMatrix" component of the discretized boundary condition has nonzero contributions if the NeumannValue has a dependency on a dependent variable; otherwise, it is zero. In other words, if the Neumann value is a generalized Neumann value, then there is a contribution to the "StiffnessMatrix" component.

The DirichletCondition representation is given through a matrix that contains positions. These positions specify which degrees of freedom should have the condition applied. A second matrix is returned that contains the values that are to be applied.

### Solution Stage

##### DeployBoundaryConditions

To have the discretized boundary conditions take effect, they need to be deployed into the system matrices. To save memory, this is done in place; no new matrices are generated.

In[53]:= |

DeployBoundaryConditions has attribute HoldFirst, so that the system matrices may be modified in place without making copies first. This means, however, that the input to DeployBoundaryConditions needs to be given in terms of symbols that have matrices as values rather than actual matrices. If the unmodified system matrices are needed later, they can be re-extracted from the DiscretizedPDEData. In other words, if the unmodified system matrices are no longer needed and computer memory consumption is crucial, then the DiscretizedPDEData needs to be cleared, e.g. with ClearAll.

During the deployment of the boundary conditions, the contributions from NeumannValue are added to the stiffness and load matrix. Then the system matrices are modified such that the DirichletCondition is satisfied.

Note that generally only after the deployment of the boundary conditions is the system of equations solvable.

##### Equation Solving

The system of equations can then be solved with LinearSolve.

### Post Processing

An interpolation function based on an ElementMesh can also be created. In this case, the second argument to NDSolve`FEM`ElementMeshInterpolation is the values at the mesh nodes.

Typically, an interpolation function extrapolates if queried outside of its domain. This behavior can be changed.

The "DefaultExtrapolationHandler" system option can be set to change the default behavior.

It is also possible to extract the ElementMesh from the interpolating function.

## Transient PDEs

Transient PDEs are time-dependent PDEs. In the simplest case, a transient PDE has a time derivative of the dependent variable as a part of the PDE. In this case, an additional system matrix is created that represents the discretized time derivative. The time integration of the system matrices can then be performed, for example, by NDSolve. Other than that the workflow is much the same as for a stationary PDE. More involved examples have coefficients and/or boundary conditions that are dependent. A later section will show how to time integrate those.

#### Transient PDE with Stationary Coefficients and Stationary Boundary Conditions—Introduction

The first example is a heat equation in 1D. Consider the following model PDE

in the spatial region 0 to 1 and a time domain from 0 to 1. Boundary and initial conditions are 0 everywhere.

In this case, the method option specifies that the spatial variable should be discretized with the finite element method and that time integration is done with the method of lines. This approach is called semi-discretization. In contrast, the full discretization would discretize both the temporal variable and the spatial variables.

The process of finding the solution to the PDE is similar to the stationary case. The main reason is that neither the PDE coefficients nor the boundary conditions depend on time. To find the solution of the transient system, the approach of semi-discretizing the PDE is taken. The PDE is discretized in space () with the finite element method. The result of that FEM discretization is then a system of ordinary differential equations in the temporal variable (). NDSolve will be used to time integrate that system of ordinary differential equations.

First, the variable data is created and populated.

Next, the solution data is set up with an initial time and a mesh.

The PDE coefficients, boundary conditions, and method data are initialized.

The discretization is performed next. Even though a transient system is considered, the discretization has only stationary components since neither the PDE coefficients nor the boundary conditions have a dependency on time.

In[74]:= |

In[75]:= |

Next, the initial conditions are set up and NDSolve will time integrate the discretized system of equations.

In[76]:= |

The NDSolve options are chosen such that they match the automatically selected options from the call to NDSolve above where NDSolve controlled all aspects of the solution process.

In[79]:= |

### Transient PDE with Stationary Coefficients and Stationary Boundary Conditions

Consider the following model PDE

subject to the boundary conditions at and at . For more information about markers, please consult the "Element Mesh Generation" tutorial.

The PDE has a first-order time derivative, but other than that the PDE coefficient and the boundary conditions are not time dependent. One strategy to solve the PDE is to first discretize the PDE and boundary conditions and then the resulting set of ODEs can be time integrated with NDSolve as a time integrator. This process is how NDSolve internally solves transient PDEs and will be shown in the next section.

For this example, a boundary mesh with markers has been predefined. The geometry models a Scanning Electrochemical Microscopy (SECM) probe. This specific probe is a heptode as it has seven electrodes embedded in a circular glass structure.

In[82]:= |

The predefined boundary mesh region has integer markers set to the boundary faces and points. Those markers can be visualized and used as positions where boundary conditions are applied during the solution process.

Since the boundary conditions only involve Dirichlet conditions, only the point elements and markers are of interest. Dirichlet conditions always operate on point elements. For completeness, the boundary mesh also has markers for the boundary elements implemented. These were needed to specify Neumann type boundary values, which always operate on boundary elements.

Strictly speaking, the node reordering is optional. When the node reordering method option is set to True, the bandwidth the system matrices have is minimized. The lower bandwidth can be beneficial for iterative solvers. Since the node reordering need only be done once, and the solvers are used for each time step, it is often advantageous to use node reordering for transient problems.

Populate the solution data. The "Time" component is the initial time t_{0} where the integration starts. Setting the "DependentVariables" and "TemporalDerivatives" components sets the initial condition and derivative of the initial condition.

In[95]:= |

Note that the first transient matrix, the damping matrix, is populated. This will be needed during the time integration later.

Note the structure of the stiffness matrix. Using "NodeReordering"->False during the mesh generation process will result in much less structured system matrices.

Before computing the transient solution, inspect the stationary solution.

In[97]:= |

In[98]:= |

In[101]:= |

Out[101]= | Play Animation |

The PDE discretization resulted in a system of ODEs. The time integration of the ODEs can be done with NDSolve. For this end, initial conditions that match the boundary conditions need to be set up.

In[103]:= |

To speed up the time integration of NDSolve, it is possible to specify the sparsity pattern of the Jacobian through an option. This is not necessary when the system is solved as a whole though NDSolve. It is necessary here because NDSolve does not try to decipher the structure for vector variables. The sparsity pattern of the Jacobian is available through the pattern of the damping matrix.

Note that using an EvaluationMonitor slows the solution process down a little bit but allows live monitoring of the solution process.

The solution shows how the transient system reaches the steady state computed above.

In[109]:= |

A single interpolating function for time and the spatial variables can be constructed. For this the time steps at which the transient interpolation has been set up need to be extracted.

Next, the interpolation values are set up. For each time step, a list of values at the nodes is required.

In[112]:= |

Next, the actual interpolation function is constructed.

### Model Order Reduction of Transient PDEs with Stationary Coefficients and Stationary Boundary Conditions

As a scenario for the transient case where access to the system matrices is important, we consider model order reduction. The question model order reduction is addressing is the following: considering a transient system of equations, is there an equivalent system of equations that is smaller than the original while preserving most of the original properties?

Consider a discretized system of equations with degrees of freedom of the form

The damping matrix and the stiffness matrix are of size and the load matrix is of size . We will now establish a second, much smaller system of equations with degrees of freedom:

Because the degrees of freedom of the second, reduced, system are much smaller than the original system , the time integration will be much more efficient. The following algorithm establishes the projection vector that maps between and :

Here is of size , is of size , and is an approximation discrepancy.

As a rough explanation, the projection vector is found as follows: Given a nonzero starting vector and a matrix , the algorithm finds orthonormal such that

In[115]:= |

The procedure returns the projection vector that has the dimensions of the reduced system times the original system.

### Transient PDEs with Transient Coefficients

The next step is to consider transient PDEs, where one or more of the coefficients is time dependent. Consider the following model PDE with a time-dependent load ,

subject to the same stationary boundary conditions as in the previous example. In this example, use a load function of

Also, variable and solution data are the same as in the previous example. The FEM method data can also be reused.

The PDE Initialization now has an additional component, the load coefficient, which it sets to the time variable .

The initialized coefficients store the coefficients specified.

In order to perform efficient time integration, the coefficients are classified into stationary and transient coefficients. Stationary coefficients are evaluated only once, much like in the previous example. Those coefficients that have time-dependent components need to be reevaluated from time step to time step and added to the stationary discretization.

The Initialization of the PDE coefficients internally performs a classification of all coefficients. This classification can be inspected and made use of. Each classification, e.g. the transient coefficients, is represented by two lists containing positions. The first list contains scalar coefficients while the second contains vector and matrix coefficients. To illustrate further, consider the general PDE form

The coefficients , , , and are scalars; , , and are vectors; and is an × matrix. The first of the two lists returned represents the scalar coefficients (, , , and ) and the second represents the vector and matrix coefficients (, , , and ).

There is a single transient component present in the model PDE.

There are no transient matrix components in this example.

The "Stationary" coefficient represents all coefficients that are stationary, including those that are not specified and hence are zero.

This classification into stationary and transient components allows for a separate computation of the stationary and transient discretization of the PDE in an efficient manner.

In the above example, the string argument "Stationary" is optional, as this is the default. Note that the load vector of the discretized stationary PDE coefficients is empty; the load vector has only transient components.

The load vector of the discretized transient PDE coefficients is populated. The time has been substituted with the numerical value 3. Also note that no other system matrices are populated, as they do not have transient coefficients.

The rationale behind the classification of PDE coefficients is that stationary coefficients can be discretized once, stored, and then reused in the time integration loop. This means that only coefficients that exhibit a time variable dependence need to be discretized during each step. This saves a lot of computational effort. Also note that since important geometric information has been precomputed and stored in the FEMMethodData object, the computational time for each transient component is kept to a minimum.

In[137]:= |

The computation of the right-hand side makes use of the stationary coefficients of the PDE that have been computed already. In this example, the damping matrix has no transient components and the factorization can be cached. If the damping matrix has time-dependent components, the above computation of the right-hand side needs to be modified, since then during each time step the damping matrix is updated and the LinearSolve step can not be cached; the factorization then needs to be performed for each time step.

Since there are no transient boundary conditions, the initial conditions are set up in the same manner as before.

In[138]:= |

The sparsity pattern is set again derived from the damping matrix. In this case, the boundary conditions have not been deployed into the damping matrix.

The visualization is done in the same manner as before.

## Coupled PDEs

This section explains how coupled partial differential equations are solved. In principle, the workflow is much the same as for single PDEs. The differences will be pointed out below.

### Deformation of a Beam under Load

As an example, a structural mechanics problem, the deformation of a beam, is considered. The deformation can be modeled as a coupled PDE in two dependent variables.

In[150]:= |

In the plane stress operator, Y is Young's modulus and ν is Poisson's ratio. Those are material properties.

The beam to be deformed is of 5 units in length and 1 unit in height.

The boundary conditions are as follows: at the left-hand side, the beam is held fixed in both x and y direction. On the right-hand boundary, a downward force of one unit is applied; a NeumannValue enforces a boundary load on the beam in the negative y direction.

In[156]:= |

The initialization has the same components as in previous sections. First, the variable and solution data are set up.

In[159]:= |

For coupled PDEs, more than one dependent variable needs to be given.

In[160]:= |

The initialization of the method data is done in the same manner as for a single dependent variable.

The display form of FEMMethodData shows the degrees of freedom, the interpolation order for each dependent variable, and integration order.

The degrees of freedom of a system of equations result from a combination of the nodes in the element mesh, the number of dependent variables, and the interpolation order of the dependent variables.

In the coupled PDE case, the coefficients become a little more involved. A coupled PDE is in essence a conglomeration of single PDEs; however, each component is present for each dependent variable. Consider the following partial differential equation in two dependent variables , :

The coupled PDE is defined in . The coefficients and are scalars in each PDE. , , and are vectors of length and is an × matrix in each PDE. The input coefficients to InitializePDECoefficients are found by correspondence of the above coupled PDE with the coupled PDE to be modeled.

In[164]:= |

In the coupled case, there is a diffusion coefficient matrix for each dependent variable in each PDE.

The display form of PDECoefficientData shows the system size, i.e. the number of dependent variables (two) and the operator size, i.e the space dimension (two) of the initialized coefficients.

Next, the boundary conditions are set up. Each of the dependent variables has a list of boundary conditions. This associates the boundary conditions with specific parts of the system of PDEs. The first set of boundary conditions is associated with the first PDE of the coupled system, the second set of boundary conditions is associated with the second set PDE of the coupled system, and so on. This allows you to uniquely associate boundary conditions with parts of PDEs.

The initialized PDE coefficients are transformed into a discrete form.

The display form of the DiscretizedPDEData shows the first dimension of the system matrices.

The visualization of the stiffness matrix shows that the structure of the underlying PDE coefficients is a system of two dependent variables. The degrees of freedom are split between the two dependent variables.

In[172]:= |

The system of equations can then be solved with LinearSolve.

Interpolation functions for each dependent variable can be created. Since the solution form LinearSolve is one solution vector containing both solutions to and , the solution vector needs to be split up.

### A Swinging Beam—Transient Coupled PDEs

For a transient coupled PDE, not much more is needed. The following example shows how to deflect the beam and have a Rayleigh damping dampen the motion out. In Rayleigh damping, the damping matrix is a linear combination of the mass matrix and the stiffness matrix. The structure of the system matrices for the two dependent variable coupled transient PDEs are represented as follows:

where the damping system matrix is a combination of the mass matrix and the stiffness matrix .

The transient coefficients also have the mass matrix coefficients set.

The display form of the DiscretizedPDEData shows the first dimension of the system matrices.

Note that the damping matrix sparse array is empty. This is because the mass matrix coefficient has been set and the default damping matrix coefficient used here is zero.

The mass matrix shows the structure of the PDE coefficients where the block diagonals are populated and the off diagonal blocks are zero.

The Rayleigh damping is computed as a linear combination of the mass and stiffness matrices.

In[185]:= |

In[186]:= |

The initial condition and the derivative of the initial condition are set to zero. This is reasonable since the Dirichlet boundary conditions are all zero and the beam starts from rest.

In[187]:= |

In[193]:= |

## Large-Scale FEM Analysis

The finite element method, as implemented in NDSolve, has been optimized for speed. This optimization comes at a cost. It is common in numerics to be able to "trade" speed versus memory consumption. This means that a faster implementation will need more memory to solve a problem at hand. It is, however, possible to tune the implemented finite element method through various means to use less memory, at the cost of some additional time that is needed to find the solution to the PDE. Some suggestions have been presented in "Solving Partial Differential Equations with Finite Elements". The following sections discuss a range of additional options to decrease the amount of memory needed during a finite element solution. The suggestions range from general advice to specific options.

For this section to be illustrative, it is good to start with a fresh kernel.

Generally, the Wolfram Language stores references to previously computed results. This mechanism can result in an accumulation of data that, if given some thought, is not needed for future computations. The mechanism that controls the history length can be set to not remember previous output.

In[195]:= |

The general procedure in the finite element discretization is as follows:

- Precompute geometric data (InitializePDEMethodData).

- Discretize the PDE (DiscretizePDE).

- Solve the system of equations (LinearSolve).

During the finite element analysis, there are two key memory bottlenecks. The first is the creation of the stiffness matrix (DiscretizePDE), and the second is the solution (LinearSolve) of the system of equations. Options to both mesh generation and geometric data precomputation have an effect on the memory requirement during discretization and solving. Furthermore, by default, the discretization step itself is performed in two steps. First all elements for a specific system matrix are computed, and then the matrix is assembled from those elements. An additional distinction is to be made if a given memory-saving action taken will merely need more time for the computation, or if that specific action also has an effect on the quality of the solution.

The options, which can also be passed to NDSolve directly, are best discussed with an example. Typical examples that require a great deal of memory are 3D application examples and/or involve many coupled PDEs.

In[196]:= |

In[197]:= |

In[198]:= |

The first crucial step is the element mesh generation. There are various aspects that influence the memory needed to solve the PDE over the mesh. The "MeshOrder" has a large effect on how much memory is needed. For example, a first-order tetrahedron has four nodes, whereas a second-order tetrahedron has 10 nodes. Reducing the "MeshOrder" also reduces the quality of the solution that can be expected, as has been discussed in "Solving Partial Differential Equations with Finite Elements".

During the matrix assembly in the upcoming discretization stage, all elements from a mesh are first computed and then assembled. This means that at one point in time the elements and the system matrix are in memory. This can be avoided by setting the option "MeshElementBlocks" to a positive integer. Then the mesh elements are split into a number of blocks. Once an element mesh block is computed, it is then assembled into the system matrix. Then the block is discarded and the next block is assembled into the same system matrix, until all blocks have been computed and assembled. With this scheme, the maximum memory consumption is the fully assembled system matrix and the size of the largest element mesh block and its computed elements. Setting "MeshElementBlocks" too large will slow down the assembly process.

During equation processing, a call to InitializePDEMethodData is made. During this call, geometric data that is needed for the finite element method is precomputed. If NDSolve`ProcessEquations is used solely for equation processing, the geometric preprocessing can be switched off.

In[208]:= |

In[209]:= |

During the discretization step, the precomputed geometry data is needed to compute the contribution of each element to the system matrices. Note that for transient systems, switching off the geometric precomputation will require that this data is recomputed in each time step. This will further increase the time needed to solve the system, but for very large systems this may be an option to obtain a solution at all.

In[217]:= |

Next, there is a very subtle point to consider. The system matrices are stored in the DiscretizedPDEData data object. Those have been extracted above and stored into the symbols load and stiffness. The deployment of the boundary conditions will modify these system matrices. But since they are stored in the DiscretizedPDEData (it is possible to extract the unmodified system matrices once more), a duplicate version will be created. To avoid this and in order to save memory, the DiscretizedPDEData object can be deleted such that only one reference to the system matrices exists.

In[221]:= |

In[224]:= |