diff --git a/utils/converters/matlab/gmsh2pdetoolbox.m b/utils/converters/matlab/gmsh2pdetoolbox.m
new file mode 100644
index 0000000000000000000000000000000000000000..d920fe2494cbe439e179fab27e16d214783daa15
--- /dev/null
+++ b/utils/converters/matlab/gmsh2pdetoolbox.m
@@ -0,0 +1,287 @@
+%GMSH2PDETOOLBOX Reads a mesh in msh format, version 1 or 2 and returns the
+%arrays p, e and t according to the MATLAB pde toolbox syntax. Filename
+%refers to the .msh file. Note that only triangular 2D mesh are processed
+%since these are the only elements allowed in pde toolbox.
+%
+%SEE ALSO load_gmsh load_gmsh4 initmesh
+%
+% Copyright (C) 2014 Arthur Levy. https://github.com/arthurlevy/
+%
+% This fucntion uses the routine load_gmsh4 in load_gmsh2.m: Copyright (C)
+% 2007 JP Moitinho de Almeida (moitinho@civil.ist.utl.pt) and R
+% Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) It is based on
+% load_gmsh.m supplied with gmsh-2.0
+function [p, e, t] = gmsh2pdetoolbox(filename)
+
+%loads the gmsh file using JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
+% and R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be) code.
+mesh_structure = load_gmsh4(filename,-1);
+
+
+%% NODES
+%nodes positions
+disp('importing nodes');
+p = mesh_structure.POS(:,1:2)';
+
+
+%% EDGES
+disp('importing edges');
+%find the edges (ie dimension 1)
+is_edge = (mesh_structure.ELE_INFOS(:,2)==1);
+%add edge data
+e = mesh_structure.ELE_NODES(is_edge,1:2)';
+e(3,:) = 0;
+e(4,:) = 1;
+%tag is important for applying boundary conditions
+e(5,:) = mesh_structure.ELE_TAGS(is_edge,1)';
+e(6,:) = 1;
+e(7,:) = 0;
+
+%% TRIANGLES
+disp('importing triangles')
+is_triangle = (mesh_structure.ELE_INFOS(:,2)==2);
+t = mesh_structure.ELE_NODES(is_triangle,1:3)';
+t(4,:) = 2;
+
+
+
+
+
+function msh = load_gmsh4(filename, which)
+
+%% Reads a mesh in msh format, version 1 or 2
+%
+% Usage:
+% To define all variables m.LINES, M.TRIANGLES, etc
+% (Warning: This creates a very large structure. Do you really need it?)
+% m = load_gmsh4('a.msh')
+%
+% To define only certain variables (for example TETS and HEXS)
+% m = load_gmsh4('a.msh', [ 5 6])
+%
+% To define no variables (i.e. if you prefer to use m.ELE_INFOS(i,2))
+% m = load_gmsh4('a.msh', -1)
+% m = load_gmsh4('a.msh', [])
+%
+% Copyright (C) 2007 JP Moitinho de Almeida (moitinho@civil.ist.utl.pt)
+% and R Lorphevre(r(point)lorphevre(at)ulg(point)ac(point)be)
+%
+% based on load_gmsh.m supplied with gmsh-2.0
+%
+% Structure msh always has the following elements:
+%
+% msh.MIN, msh.MAX - Bounding box
+% msh.nbNod - Number of nodes
+% msh.nbElm - Total number of elements
+% msh.nbType(i) - Number of elements of type i (i in 1:19)
+% msh.POS(i,j) - j'th coordinate of node i (j in 1:3, i in 1:msh.nbNod)
+% msh.ELE_INFOS(i,1) - id of element i (i in 1:msh.nbElm)
+% msh.ELE_INFOS(i,2) - type of element i
+% msh.ELE_INFOS(i,3) - number of tags for element i
+% msh.ELE_INFOS(i,4) - Dimension (0D, 1D, 2D or 3D) of element i
+% msh.ELE_TAGS(i,j) - Tags of element i (j in 1:msh.ELE_INFOS(i,3))
+% msh.ELE_NODES(i,j) - Nodes of element i (j in 1:k, with
+% k = msh.NODES_PER_TYPE_OF_ELEMENT(msh.ELE_INFOS(i,2)))
+%
+% These elements are created if requested:
+%
+% msh.nbLines - number of 2 node lines
+% msh.LINES(i,1:2) - nodes of line i
+% msh.LINES(i,3) - tag (WHICH ?????) of line i
+%
+% msh.nbTriangles - number of 2 node triangles
+% msh.TRIANGLES(i,1:3) - nodes of triangle i
+% msh.TRIANGLES(i,4) - tag (WHICH ?????) of triangle i
+%
+% Etc
+
+% These definitions need to be updated when new elemens types are added to gmsh
+%
+% msh.Types{i}{1} Number of an element of type i
+% msh.Types{i}{2} Dimension (0D, 1D, 2D or 3D) of element of type i
+% msh.Types{i}{3} Name to add to the structure with elements of type i
+% msh.Types{i}{4} Name to add to the structure with the number of elements of type i
+%
+
+nargchk(1, 2, nargin);
+
+msh.Types = { ...
+ { 2, 1, 'LINES', 'nbLines'}, ... % 1
+ { 3, 2, 'TRIANGLES', 'nbTriangles'}, ...
+ { 4, 2, 'QUADS', 'nbQuads'}, ...
+ { 4, 3, 'TETS', 'nbTets'}, ...
+ { 8, 3, 'HEXAS', 'nbHexas'}, ... %5
+ { 6, 3, 'PRISMS', 'nbPrisms'}, ...
+ { 5, 3, 'PYRAMIDS', 'nbPyramids'}, ...
+ { 3, 1, 'LINES3', 'nbLines3'}, ...
+ { 6, 2, 'TRIANGLES6', 'nbTriangles6'}, ...
+ { 9, 2, 'QUADS9', 'nbQuads9'}, ... % 10
+ { 10, 3, 'TETS10', 'nbTets10'}, ...
+ { 27, 3, 'HEXAS27', 'nbHexas27'}, ...
+ { 18, 3, 'PRISMS18', 'nbPrisms18'}, ...
+ { 14, 3, 'PYRAMIDS14', 'nbPyramids14'}, ...
+ { 1, 0, 'POINTS', 'nbPoints'}, ... % 15
+ { 8, 3, 'QUADS8', 'nbQuads8'}, ...
+ { 20, 3, 'HEXAS20', 'nbHexas20'}, ...
+ { 15, 3, 'PRISMS15', 'nbPrisms15'}, ...
+ { 13, 3, 'PYRAMIDS13', 'nbPyramids13'}, ...
+};
+
+ntypes = length(msh.Types);
+
+if (nargin==1)
+ which = 1:ntypes;
+else
+ if isscalar(which) && which == -1
+ which = [];
+ end
+end
+
+% Could check that "which" is properlly defined....
+
+fid = fopen(filename, 'r');
+fileformat = 0; % Assume wrong file
+
+tline = fgetl(fid);
+if (feof(fid))
+ disp (sprintf('Empty file: %s', filename));
+ msh = [];
+ return;
+end
+
+%% Read mesh type
+if (strcmp(tline, '$NOD'))
+ fileformat = 1;
+elseif (strcmp(tline, '$MeshFormat'))
+ fileformat = 2;
+ tline = fgetl(fid);
+ if (feof(fid))
+ disp (sprintf('Syntax error (no version) in: %s', filename));
+ fileformat = 0;
+ end
+ [ form ] = sscanf( tline, '%f %d %d');
+% if (form(1) ~= 2)
+% disp (sprintf('Unknown mesh format: %s', tline));
+% fileformat = 0;
+% end
+ if (form(2) ~= 0)
+ disp ('Sorry, this program can only read ASCII format');
+ fileformat = 0;
+ end
+ fgetl(fid); % this should be $EndMeshFormat
+ if (feof(fid))
+ disp (sprintf('Syntax error (no $EndMeshFormat) in: %s', filename));
+ fileformat = 0;
+ end
+ tline = fgetl(fid); % this should be $Nodes
+ if (feof(fid))
+ disp (sprintf('Syntax error (no $Nodes) in: %s', filename));
+ fileformat = 0;
+ end
+end
+
+if (~fileformat)
+ msh = [];
+ return
+end
+
+%% Read nodes
+
+if strcmp(tline, '$NOD') || strcmp(tline, '$Nodes')
+ msh.nbNod = fscanf(fid, '%d', 1);
+ aux = fscanf(fid, '%g', [4 msh.nbNod]);
+ msh.POS = aux(2:4,:)';
+ numids = max(aux(1,:));
+ IDS = zeros(1, numids);
+ IDS(aux(1,:)) = 1:msh.nbNod; % This may not be an identity
+ msh.MAX = max(msh.POS);
+ msh.MIN = min(msh.POS);
+ fgetl(fid); % End previous line
+ fgetl(fid); % Has to be $ENDNOD $EndNodes
+else
+ disp (sprintf('Syntax error (no $Nodes/$NOD) in: %s', filename));
+ fileformat = 0;
+end
+
+%% Read elements
+tline = fgetl(fid);
+if strcmp(tline,'$ELM') || strcmp(tline, '$Elements')
+ msh.nbElm = fscanf(fid, '%d', 1);
+ % read all info about elements into aux (it is faster!)
+ aux = fscanf(fid, '%g', inf);
+ start = 1;
+ msh.ELE_INFOS = zeros(msh.nbElm, 4); % 1 - id, 2 - type, 3 - ntags, 4 - Dimension
+ msh.ELE_NODES = zeros(msh.nbElm,6); % i - Element, j - ElNodes
+ if (fileformat == 1)
+ ntags = 2;
+ else
+ ntags = 3; % just a prediction
+ end
+ msh.ELE_TAGS = zeros(msh.nbElm, ntags); % format 1: 1 - physical number, 2 - geometrical number
+ % format 2: 1 - physical number, 2 - geometrical number, 3 - mesh partition number
+ msh.nbType = zeros(ntypes,1);
+ for I = 1:msh.nbElm
+ if (fileformat == 2)
+ finnish = start + 2;
+ msh.ELE_INFOS(I, 1:3) = aux(start:finnish);
+ ntags = aux(finnish);
+ start = finnish + 1;
+ finnish = start + ntags -1;
+ msh.ELE_TAGS(I, 1:ntags) = aux(start:finnish);
+ start = finnish + 1;
+ else
+ finnish = start + 1;
+ msh.ELE_INFOS(I, 1:2) = aux(start:finnish);
+ start = finnish + 1; % the third element is nnodes, which we get from the type
+ msh.ELE_INFOS(I, 3) = 2;
+ finnish = start + 1;
+ msh.ELE_TAGS(I, 1:2) = aux(start:finnish);
+ start = finnish + 2;
+ end
+ type = msh.ELE_INFOS(I, 2);
+ msh.nbType(type) = msh.nbType(type) + 1;
+ msh.ELE_INFOS(I, 4) = msh.Types{type}{2};
+ nnodes = msh.Types{type}{1};
+ finnish = start + nnodes - 1;
+ msh.ELE_NODES(I, 1:nnodes) = IDS(aux(start:finnish));
+ start=finnish + 1;
+ end
+ fgetl(fid); % Has to be $ENDELM or $EndElements
+else
+ disp (sprintf('Syntax error (no $Elements/$ELM) in: %s', filename));
+ fileformat = 0;
+end
+
+if (fileformat == 0)
+ msh = [];
+end
+
+fclose(fid);
+
+%% This is used to create the explicit lists for types of elements
+
+for i = which
+ if (~isempty(msh.Types{i}{3}))
+ cmd = sprintf('msh.%s=msh.nbType(%d);', msh.Types{i}{4}, i);
+ eval(cmd);
+ % Dimension
+ cmd = sprintf('msh.%s=zeros(%d,%d);', msh.Types{i}{3}, msh.nbType(i), msh.Types{i}{1}+1);
+ eval(cmd);
+ % Clear nbType for counting, next loop will recompute it
+ msh.nbType(i) = 0;
+ end
+end
+
+for i = 1:msh.nbElm
+ type = msh.ELE_INFOS(i,2);
+ if (find(which == type))
+ if (~isempty(msh.Types{type}{3}))
+ msh.nbType(type) = msh.nbType(type)+1;
+ aux=[ msh.ELE_NODES(i,1:msh.Types{type}{1}), msh.ELE_TAGS(i,1) ];
+ cmd = sprintf('msh.%s(%d,:)=aux;', msh.Types{type}{3}, msh.nbType(type));
+ eval(cmd);
+ end
+ end
+end
+
+return;