Utilities

MOLE provides a set of utility functions and classes to simplify common tasks when working with mimetic operators and boundary conditions.

API Reference

class Utils

Utility Functions.

Public Functions

void meshgrid(const vec &x, const vec &y, mat &X, mat &Y)

An analog to the MATLAB/Octave 2D meshgrid operation.

returns 2-D grid coordinates based on the coordinates contained in vectors x and y. X is a matrix where each row is a copy of x, and Y is a matrix where each column is a copy of y. The grid represented by the coordinates X and Y has length(y) rows and length(x) columns. Key here is the rows is the y-coordinate, and the columns are the x-coordinate.

Parameters:

  • x – a vector of x-indices

  • y – a vector of y-indices

  • X – a sparse matrix, will be filled by the function

  • Y – a sparse matrix, will be filled by the function

void meshgrid(const vec &x, const vec &y, const vec &z, cube &X, cube &Y, cube &Z)

An analog to the MATLAB/Octave 3D meshgrid operation.

meshgrid(x,y,z,X,Y,Z) returns 3-D grid coordinates defined by the vectors x, y, and z. The grid represented by X, Y, and Z has size length(y)-by-length(x)-by-length(z).

Parameters:

  • x – a vector of x-indices

  • y – a vector of y-indices

  • z – a vector of z-indices

  • X – a sparse matrix, will be filled by the function

  • Y – a sparse matrix, will be filled by the function

  • Z – a sparse matrix, will be filled by the function

Public Static Functions

static sp_mat spkron(const sp_mat &A, const sp_mat &B)

A wrappper for implementing a sparse Kroenecker product.

Note

This is available in Armadillo >8.0

Parameters:

  • A – a sparse matrix

  • B – a sparse matrix

static sp_mat spjoin_rows(const sp_mat &A, const sp_mat &B)

An in place operation for joining two matrices by rows.

Note

This is available in Armadillo >8.0

Parameters:

  • A – a sparse matrix

  • B – a sparse matrix

static sp_mat spjoin_cols(const sp_mat &A, const sp_mat &B)

An in place operation for joining two matrices by columns.

Note

This is available in Armadillo >=8.5

Parameters:

  • A – a sparse matrix

  • B – a sparse matrix

static vec spsolve_eigen(const sp_mat &A, const vec &b)

A wrappper for implementing a sparse solve using Eigen from SuperLU.

Note

This function requires the EIGEN to be used when Armadillo is built

Parameters:

  • A – a sparse matrix LHS of Ax=b

  • b – a vector for the RHS of Ax=b

static double trapz(const vec &x, const vec &y)

Implements the trapezoidal rule for 1D numerical integration.

Computes the area under a curve defined by vectors x and y using: A ≈ ∑ 0.5 * (xᵢ₊₁ - xᵢ) * (yᵢ + yᵢ₊₁)

Parameters:

  • x – Vector of x-coordinates

  • y – Vector of y-values at corresponding x

Returns:

Estimated area under the curve

Usage Examples

Here’s an example using utility functions in a parabolic equation:

Parabolic 1D Example (examples/cpp/parabolic1D.cpp)
 1#include <iostream>
 2#include <math.h> 
 3#include "mole.h"
 4/**
 5 * This example uses MOLE to solve the heat equation u_t-alpha*u_xx=0[with alpha=1] over [0,1]^2,
 6 * u(x,0)=0 for x in (-1,1), u(1,t)=u(-1,t)=100 for t in [0,1]. 
 7 */
 8int main() { 
 9    int k=2;// Operators' order of accuracy
10    double t0=0; //initial time
11    double tf=1; //final time
12    double a=0; //left boundary
13    double b=1; //right boundary
14    int m=2*k+1; //num of cells
15    double dx=(b-a)/m;
16    double dt=tf/(ceil((3*tf)/(dx*dx))); //Von Neumann stability criterion for explicit scheme, if k > 2 then dt/dx^2<0.5.  
17    /*Note that unlike the matlab example, dt isn't dx^2/3 because if a and b are changed or the final time is changed
18    *the final time tf may not be a multiple of dt. This value of dt guarantees that the final time will be a mutiple of dt.
19    * while still satisfying Von Neumann stablity criterion
20    *Note that in this example dt=dx^2/3 like it is in the matlab example. 
21    */
22    Laplacian L(k,m,dx); 
23    vec solution(m + 2); 
24    solution(0)=100;
25    solution(m+1)=100;
26    vec k1(m+2); 
27    double t=t0;
28    while (t <= tf) { //time integration with euler method. 
29        k1=L*(solution);
30        solution=solution+dt*k1;
31        t=t+dt;
32    }
33    std::cout << solution;
34
35
36    return 0;
37}