From The FEniCS project
The power and efficiency of DOLFIN can be integrated into the familiar and flexible Matlab environment. One smart and simple way to do this involves compiling functions as Matlab EXecutables (MEX). Compiling a C/Fortran program with the mex-compiler results in an object that is called as an ordinary Matlab function. The objective here is to show how to wrap a Poisson solver (from the DOLFIN demos) so that it can be called from within Matlab, like
> [u vert tri edge] = poisson_DOLFIN('dolfin-2.xml.gz')
> pdeplot(vert, edg, tri,'xydata',U,'zdata',U,'mesh','on')
Contents |
MEX Basics
The C-MEX API is well documented by The MathWorks:
Any MEX program provides an entry-point, namely mexFunction, as illustrated by this "Hello World" program:
#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[])
{
cout << "Hello, world! << endl;
}
Some notes about the components here:
- int nrhs is the number of input arguments
- int nrls is the number of output arguments
- mxArray is the basic (and only) data structure that provides the interface for numeric data input from (prhs[]), and output to (plhs[]), Matlab.
MEX file for Poisson solver
Download the tarball that contains this example here
The main job of the MEX program is to allocate data for the results of the computation, and then call the Poisson solver (which is completely separated from the MEX-specifics).
In its most basic form, we could write something like
#include "mex.h"
#include "psolver.h"
// Input args
#define N prhs[0] // N-by-N UnitSquare mesh
// Output destinations
#define U plhs[0] // Solution vector
using namespace dolfin;
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{
// Initialize DOLFIN-based solver
PoissonSolver s( (int)*mxGetPr(N) );
// Get mesh size and allocate data
int num_vert = s.getNumVertices();
U = mxCreateDoubleMatrix(num_vert, 1, mxREAL);
// Compute solution
s.solve(mxGetPr(U));
}
However, we want to return the mesh also, and have as input either a mesh file or a integer that specifies the size of a UnitSquare mesh. The following code is found in the file poisson_DOLFIN.cpp:
#include "mex.h"
#include "psolver.h"
// Input args
#define GRID prhs[0] // mesh file or NxN unit square
// Output destinations
#define U plhs[0] // Solution vector
#define VERTEX plhs[1] // Mesh vertices
#define TRI plhs[2] // Mesh triangels
#define EDGE plhs[3] // Mesh edges (dummy)
using namespace dolfin;
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{
// Parse input
if (nrhs != 1)
mexErrMsgTxt("Poisson_mex: Must have one input argument. ");
else if (nlhs =! 4)
mexErrMsgTxt("Poisson_mex: Must have four output arguments.");
// Initialize DOLFIN-based solver
PoissonSolver *s;
if (mxIsChar(GRID))
{
int buflen = mxGetNumberOfElements(GRID) + 1;
char* buf = (char*) mxCalloc(buflen, sizeof(char));
mxGetString(GRID, buf, buflen);
s = new PoissonSolver( std::string(buf) );
}
else
s = new PoissonSolver( (int)*mxGetPr(GRID) );
// Get mesh size and allocate data
int num_vert = s->getNumVertices();
U = mxCreateDoubleMatrix(num_vert, 1, mxREAL);
// Compute solution
s->solve(mxGetPr(U));
// Get mesh
int num_tri = s->getNumCells();
int num_edges = s->getNumEdges();
VERTEX = mxCreateDoubleMatrix(2, num_vert, mxREAL);
TRI = mxCreateNumericMatrix(4, num_tri, mxINT32_CLASS, 0);
EDGE = mxCreateNumericMatrix(7, num_edges, mxINT32_CLASS,0);
s->getMesh(mxGetPr(VERTEX), (int*)mxGetData(TRI), (int*)mxGetData(EDGE));
}
The PoissonSolver class is a rearrangement of the Poisson demo, with some added functionality for exporting the mesh. Note that the edge output is a dummy (since edge data is not really necessary for plotting in Matlab with pdeplot). Instead of writing data to disk, the computed solution is interpolated at vertices and stored in the output vector, i.e.
void PoissonSolver::solve(double* solution)
{
(...)
Function u;
pde.solve(u);
// Write solution to mxArray
u.interpolate(solution);
}
Compiling
Matlab must be set up to use a C++ compiler (I used gcc). This is done by calling
> mex -setup
Then the Poisson solver can be compiled and linked with
> mex poisson_DOLFIN.cpp psolver.cpp -ldolfin -lamd -lumfpack -lblas -lxml2
To see that all went well, run
> [U vert tri edg] = poisson_DOLFIN(32); > pdeplot(vert, double(edg), double(tri),'xydata',U,'zdata',U,'mesh','on')
You should get the familiar Poisson solution from the basic demo:
Potential problems
On a recent Ubuntu box you might get a run-time abort from the loader with an error message like
/usr/lib/libxml2.so.2: undefined symbol: gzopen64
To diagnose the problem, inspect the following libraries:
host:poisson > objdump -T /usr/lib/libxml2.so.2 | grep gzopen 00000000 DF *UND* 0000001b gzopen64 host:poisson > objdump -T /usr/lib/libz.so.1 | grep gzopen 00004060 g DF .text 0000001b Base gzopen 00004040 g DF .text 0000001b Base gzopen64 host:poisson > objdump -T <matlab path>/bin/<arch>/libz.so.1 | grep gzopen 00003de0 g DF .text 00000014 Base gzopen
From this we learn that libxml2 makes a reference to gzopen64, which is found in the system-wide libz but not in the libz provided by the Matlab distribution. The problem that occurs is that Matlab loads its own version of libz, which lacks the instructions required by libxml2. One solution to this problem is simply to move Matlab's internal libz* to a back-up folder.
