Index: /issm/trunk-jpl/src/m/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/m/Makefile.am	(revision 22876)
+++ /issm/trunk-jpl/src/m/Makefile.am	(revision 22877)
@@ -82,77 +82,21 @@
 if !DEVELOPMENT
 
-js_scripts = ${ISSM_DIR}/src/m/array/arrayoperations.js \
-				${ISSM_DIR}/src/m/boundaryconditions/SetIceSheetBC.js \
-				${ISSM_DIR}/src/m/boundaryconditions/SetIceShelfBC.js \
-				${ISSM_DIR}/src/m/classes/amr.js \
-				${ISSM_DIR}/src/m/classes/autodiff.js \
-				${ISSM_DIR}/src/m/classes/balancethickness.js \
-				${ISSM_DIR}/src/m/classes/basalforcings.js \
-				${ISSM_DIR}/src/m/classes/calving.js \
-				${ISSM_DIR}/src/m/classes/clusters/generic.js \
-				${ISSM_DIR}/src/m/classes/clusters/local.js \
-				${ISSM_DIR}/src/m/classes/constants.js \
-				${ISSM_DIR}/src/m/classes/damage.js \
-				${ISSM_DIR}/src/m/classes/debug.js \
-				${ISSM_DIR}/src/m/classes/flowequation.js \
-				${ISSM_DIR}/src/m/classes/friction.js \
-				${ISSM_DIR}/src/m/classes/geometry.js \
-				${ISSM_DIR}/src/m/classes/giaivins.js \
-				${ISSM_DIR}/src/m/classes/fourierlove.js \
-				${ISSM_DIR}/src/m/classes/esa.js \
-				${ISSM_DIR}/src/m/classes/groundingline.js \
-				${ISSM_DIR}/src/m/classes/hydrologyshreve.js \
-				${ISSM_DIR}/src/m/classes/initialization.js \
-				${ISSM_DIR}/src/m/classes/inversion.js \
-				${ISSM_DIR}/src/m/classes/mask.js \
-				${ISSM_DIR}/src/m/classes/maskpsl.js \
-				${ISSM_DIR}/src/m/classes/masstransport.js \
-				${ISSM_DIR}/src/m/classes/matice.js \
-				${ISSM_DIR}/src/m/classes/mesh2d.js \
-				${ISSM_DIR}/src/m/classes/mesh3dprisms.js \
-				${ISSM_DIR}/src/m/classes/mesh3dsurface.js \
-				${ISSM_DIR}/src/m/classes/miscellaneous.js \
-				${ISSM_DIR}/src/m/classes/model.js \
-				${ISSM_DIR}/src/m/classes/outputdefinition.js \
-				${ISSM_DIR}/src/m/classes/pairoptions.js \
-				${ISSM_DIR}/src/m/classes/plotoptions.js \
-				${ISSM_DIR}/src/m/classes/priv.js \
-				${ISSM_DIR}/src/m/classes/qmu.js \
-				${ISSM_DIR}/src/m/classes/radaroverlay.js \
-				${ISSM_DIR}/src/m/classes/rifts.js \
-				${ISSM_DIR}/src/m/classes/slr.js \
-				${ISSM_DIR}/src/m/classes/issmsettings.js \
-				${ISSM_DIR}/src/m/classes/SMBforcing.js \
-				${ISSM_DIR}/src/m/classes/steadystate.js \
-				${ISSM_DIR}/src/m/classes/stressbalance.js \
-				${ISSM_DIR}/src/m/classes/template.js \
-				${ISSM_DIR}/src/m/classes/thermal.js \
-				${ISSM_DIR}/src/m/classes/timestepping.js \
-				${ISSM_DIR}/src/m/classes/toolkits.js \
-				${ISSM_DIR}/src/m/classes/levelset.js \
-				${ISSM_DIR}/src/m/classes/trans.js \
-				${ISSM_DIR}/src/m/classes/verbose.js \
-				${ISSM_DIR}/src/m/consistency/checkfield.js \
-				${ISSM_DIR}/src/m/consistency/ismodelselfconsistent.js \
+js_scripts = ${ISSM_DIR}/src/m/array/*.js \
+				${ISSM_DIR}/src/m/boundaryconditions/*.js \
+				${ISSM_DIR}/src/m/classes/*.js \
+                ${ISSM_DIR}/src/m/classes/clusters/*.js \
+				${ISSM_DIR}/src/m/consistency/*.js \
 				${ISSM_DIR}/src/m/extrusion/*.js \
-				${ISSM_DIR}/src/m/geometry/FlagElements.js \
-				${ISSM_DIR}/src/m/inversions/marshallcostfunctions.js \
+				${ISSM_DIR}/src/m/geometry/*.js \
+				${ISSM_DIR}/src/m/inversions/*.js \
 				${ISSM_DIR}/src/m/io/*.js \
-				${ISSM_DIR}/src/m/materials/paterson.js \
-				${ISSM_DIR}/src/m/mesh/triangle.js \
-				${ISSM_DIR}/src/m/miscellaneous/colorbars.js \
-				${ISSM_DIR}/src/m/miscellaneous/fielddisplay.js \
-				${ISSM_DIR}/src/m/miscellaneous/rgbcolor.js \
-				${ISSM_DIR}/src/m/parameterization/setflowequation.js \
-				${ISSM_DIR}/src/m/parameterization/setmask.js \
+				${ISSM_DIR}/src/m/materials/*.js \
+				${ISSM_DIR}/src/m/mesh/*.js \
+				${ISSM_DIR}/src/m/miscellaneous/*.js \
+				${ISSM_DIR}/src/m/parameterization/*.js \
 				${ISSM_DIR}/src/m/plot/*.js \
-				${ISSM_DIR}/src/m/print/sprintf.js \
-				${ISSM_DIR}/src/m/solve/loadresultsfrombuffer.js \
-				${ISSM_DIR}/src/m/solve/marshall.js \
-				${ISSM_DIR}/src/m/solve/parseresultsfrombuffer.js \
-				${ISSM_DIR}/src/m/solve/solve.js \
-				${ISSM_DIR}/src/m/solve/WriteData.js \
-				${ISSM_DIR}/src/m/solvers/issmgslsolver.js \
-				${ISSM_DIR}/src/m/solvers/issmmumpssolver.js
+				${ISSM_DIR}/src/m/print/*.js \
+				${ISSM_DIR}/src/m/solve/*.js \
+				${ISSM_DIR}/src/m/solvers/*.js
 bin_SCRIPTS +=  issm-bin.js
 
Index: /issm/trunk-jpl/src/m/array/arrayoperations.js
===================================================================
--- /issm/trunk-jpl/src/m/array/arrayoperations.js	(revision 22876)
+++ /issm/trunk-jpl/src/m/array/arrayoperations.js	(revision 22877)
@@ -47,4 +47,100 @@
 	for(var i=0;i<array.length;i++)sum+=array[i];
 	return sum;
+} //}}}
+
+function ArrayAdd(array1, array2){ //{{{
+    if(array1.length!=array2.length)throw Error("ArrayAdd error message: arrays provided as arguments are not of the same length!");
+
+    //do the subtraction:
+    var result = array1.slice();
+    for(var i = 0; i < array1.length; i++){
+        result[i] = array1[i] + array2[i];
+    }
+    return result;
+
+} //}}}
+function ArraySubtract(array1, array2){ //{{{
+    if(array1.length!=array2.length)throw Error("ArraySubtract error message: arrays provided as arguments are not of the same length!");
+
+    //do the subtraction:
+    var result = array1.slice();
+    for(var i = 0; i < array1.length; i++){
+        result[i] = array1[i] - array2[i];
+    }
+    return result;
+
+} //}}}
+function ArrayMultiply(array1, array2){ //{{{
+    if(array1.length!=array2.length)throw Error("ArrayMultiply error message: arrays provided as arguments are not of the same length!");
+
+    //do the subtraction:
+    var result = array1.slice();
+    for(var i = 0; i < array1.length; i++){
+        result[i] = array1[i] * array2[i];
+    }
+    return result;
+
+} //}}}
+function ArrayDivide(array1, array2){ //{{{
+    if(array1.length!=array2.length)throw Error("ArrayDivide error message: arrays provided as arguments are not of the same length!");
+
+    //do the subtraction:
+    var result = array1.slice();
+    for(var i = 0; i < array1.length; i++){
+        result[i] = array1[i] / array2[i];
+    }
+    return result;
+
+} //}}}
+function ArrayMean(array){ //{{{
+    //Calculate the mean:
+    var sum = ArraySum(array);
+    return sum / array.length;
+
+} //}}}
+function ArraySTD(array){ //{{{
+    //Calculate the standard deviation:
+    var sum = ArraySum(array);
+    var differences = 0;
+    for(var i = 0; i < array.length; i++){
+        differences += (array[i] - sum) ** 2;
+    }
+    var variance = differences / array.length;
+    return Math.sqrt(variance);
+
+} //}}}
+function ArrayLerp(array1, array2, weight1, weight2, alpha){ //{{{
+    //Returns linear combination of array1 and array2, based on lerp coefficients weight1, weight2, and alpha.
+    //Parameters:
+    //  array1  - array 
+    //  array2  - array 
+    //  weight1 - value associated with array1.
+    //  weight2 - value associated with array2
+    //  alpha   - value associated with output array
+    //Output:
+    //  result - linear combination of array1 and array2 based on the equation array1 * alpha + 
+    
+    var range = weight2 - weight1;
+    var normalizedAlpha = (alpha - weight1) / range;
+    var result = array1.slice();
+    for(var i = 0; i < array1.length; i++){
+        result[i]  = array1[i] * normalizedAlpha + array2[i] * (1 - normalizedAlpha);
+    }
+    return result;
+
+} //}}}
+function ArrayTranspose(array){ //{{{
+    //Transposes 2d array.
+    var rows = array.length;
+    var cols = array[0].length;
+    var result = Create2DArray(cols, rows);
+
+    for(var i = 0; i < rows; i++) {
+        for(var j = 0; j < cols; j++) {
+            result[j][i]  = array[i][j];
+        }
+    }
+    return result;
+
 } //}}}
 function ArrayXPY(){ //{{{
@@ -93,11 +189,4 @@
 	return or;
 
-} //}}}
-function ArrayTranspose(array) { //{{{
-	return Object.keys(array[0]).map(function(column) {
-		return array.map(function(row) { 
-			return row[column]; 
-		});
-	});
 } //}}}
 function ArrayConcat(a, b) { //{{{
@@ -227,4 +316,12 @@
 	});
 
+} //}}}
+function ArrayAny(array) { //{{{
+    //Emulates Matlab 'any' function
+    
+    for(var i=0;i<array.length;i++){
+        if (array[i]!=0)return 1;
+    }
+    return 0;
 } //}}}
 function ArrayAnyEqual(array,value) { //{{{
Index: /issm/trunk-jpl/src/m/classes/bamggeom.js
===================================================================
--- /issm/trunk-jpl/src/m/classes/bamggeom.js	(revision 22877)
+++ /issm/trunk-jpl/src/m/classes/bamggeom.js	(revision 22877)
@@ -0,0 +1,28 @@
+//BAMGGEOM class definition
+//
+//   Usage:
+//      bamggeom(varargin)
+
+function bamggeom(){
+    //methods
+    this.setdefaultparameters = function(){// {{{
+    }// }}}
+    this.disp= function(){// {{{
+        disp(sprintf('\n%s = \n', 'bamggeom'));
+        disp(this);
+    }// }}}
+
+    //properties 
+    // {{{
+    this.Vertices           = [];
+    this.Edges              = [];
+    this.TangentAtEdges     = [];
+    this.Corners            = [];
+    this.RequiredVertices   = [];
+    this.RequiredEdges      = [];
+    this.CrackedEdges       = [];
+    this.SubDomains         = [];
+
+    this.setdefaultparameters();
+    //}}}
+}
Index: /issm/trunk-jpl/src/m/classes/bamgmesh.js
===================================================================
--- /issm/trunk-jpl/src/m/classes/bamgmesh.js	(revision 22877)
+++ /issm/trunk-jpl/src/m/classes/bamgmesh.js	(revision 22877)
@@ -0,0 +1,36 @@
+//BAMGMESH class definition
+//
+//   Usage:
+//      bamgmesh(varargin)
+
+function bamgmesh(){
+    //methods
+    this.setdefaultparameters = function(){// {{{
+    }// }}}
+    this.disp= function(){// {{{
+        disp(sprintf('\n%s = \n', 'bamgmesh'));
+        disp(this);
+    }// }}}
+
+    //properties 
+    // {{{
+    this.Vertices                   = [];
+    this.Edges                      = [];
+    this.Triangles                  = [];
+    this.IssmEdges                  = [];
+    this.IssmSegments               = [];
+    this.VerticesOnGeomVertex       = [];
+    this.VerticesOnGeomEdge         = [];
+    this.EdgesOnGeomEdge            = [];
+    this.SubDomains                 = [];
+    this.SubDomainsFromGeom         = [];
+    this.ElementConnectivity        = [];
+    this.NodalConnectivity          = [];
+    this.NodalElementConnectivity   = [];
+    this.CrackedVertices            = [];
+    this.CrackedEdges               = [];
+    this.PreviousNumbering          = [];
+
+    this.setdefaultparameters();
+    //}}}
+}
Index: /issm/trunk-jpl/src/m/mesh/bamg.js
===================================================================
--- /issm/trunk-jpl/src/m/mesh/bamg.js	(revision 22877)
+++ /issm/trunk-jpl/src/m/mesh/bamg.js	(revision 22877)
@@ -0,0 +1,701 @@
+function bamg(md){
+    //BAMG - mesh generation
+    //
+    //   Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):
+    //
+    //   - domain :            followed by an ARGUS file that prescribes the domain outline
+    //   - holes :             followed by an ARGUS file that prescribes the holes
+    //   - subdomains :        followed by an ARGUS file that prescribes the list of
+    //                         subdomains (that need to be inside domain)
+    //
+    //   - hmin :              minimum edge length (default is 10^-100)
+    //   - hmax :              maximum edge length (default is 10^100)
+    //   - hVertices :         imposed edge length for each vertex (geometry or mesh)
+    //   - hminVertices :      minimum edge length for each vertex (mesh)
+    //   - hmaxVertices :      maximum edge length for each vertex (mesh)
+    //
+    //   - anisomax :          maximum ratio between the smallest and largest edges (default is 10^30)
+    //   - coeff :             coefficient applied to the metric (2-> twice as many elements, default is 1)
+    //   - cutoff :            scalar used to compute the metric when metric type 2 or 3 are applied
+    //   - err :               error used to generate the metric from a field
+    //   - errg :              geometric error (default is 0.1)
+    //   - field :             field of the model that will be used to compute the metric
+    //                         to apply several fields, use one column per field
+    //   - gradation :         maximum ratio between two adjacent edges
+    //   - Hessiantype :       0 -> use double L2 projection (default)
+    //                         1 -> use Green formula
+    //   - KeepVertices :      try to keep initial vertices when adaptation is done on an existing mesh (default 1)
+    //   - NoBoundaryRefinment: do not refine boundary, only follow contour provided (default 0)
+    //   - maxnbv :            maximum number of vertices used to allocate memory (default is 10^6)
+    //   - maxsubdiv :         maximum subdivision of exisiting elements (default is 10)
+    //   - metric :            matrix (numberofnodes x 3) used as a metric
+    //   - Metrictype :        0 -> absolute error          c/(err coeff^2) * Abs(H)        (default)
+    //                         1 -> relative error          c/(err coeff^2) * Abs(H)/max(s,cutoff*max(s))
+    //                         2 -> rescaled absolute error c/(err coeff^2) * Abs(H)/(smax-smin)
+    //   - nbjacoby :          correction used by Hessiantype=1 (default is 1)
+    //   - nbsmooth :          number of metric smoothing procedure (default is 3)
+    //   - omega :             relaxation parameter of the smoothing procedure (default is 1.8)
+    //   - power :             power applied to the metric (default is 1)
+    //   - splitcorners :      split triangles whuch have 3 vertices on the outline (default is 1)
+    //   - verbose :           level of verbosity (default is 1)
+    //
+    //   - vertical :          is this a 2d vertical mesh (flowband, default is 0)
+    //   - rifts :             followed by an ARGUS file that prescribes the rifts
+    //   - toltip :            tolerance to move tip on an existing point of the domain outline
+    //   - tracks :            followed by an ARGUS file that prescribes the tracks that the mesh will stick to
+    //   - RequiredVertices :  mesh vertices that are required. [x,y,ref]; ref is optional
+    //   - tol :               if the distance between 2 points of the domain outline is less than tol, they
+    //                         will be merged
+    //
+    //   Examples:
+    //      md=bamg(md,'domain','DomainOutline.exp','hmax',3000);
+    //      md=bamg(md,'field',[md.inversion.vel_obs md.geometry.thickness],'hmax',20000,'hmin',1000);
+    //      md=bamg(md,'metric',A,'hmin',1000,'hmax',20000,'gradation',3,'anisomax',1);
+
+    //process options
+    var args = Array.prototype.slice.call(arguments);
+    var options = new pairoptions(args.slice(1,args.length));
+    options.deleteduplicates(1);
+
+    //initialize the structures required as input of Bamg
+    var bamg_options = {}
+    var bamg_geometry = bamggeom();
+    var bamg_mesh = bamgmesh();
+
+    var subdomain_ref = 1;
+    var hole_ref = 1;
+    var vertices = [];
+    var edges = [];
+    var subdomains = [];
+    // Bamg Geometry parameters {{{
+    if (options.exist('domain')) {
+
+        //Check that file exists
+        var domainfile=options.getfieldvalue('domain');
+        if ((typeof domainfile) === 'string') {
+            console.log('bamg error message: file ' + domainfile + ' loading from file path not supported - domain must be file object (.shp or .exp)');
+        } else if ((typeof domainfile) === 'object') {
+            domain = domainfile;
+        } else {
+            console.log('"domain" type not supported yet - ' + (typeof domainfile));
+        }
+
+        var holes = [];
+        if (options.exist('holes')) {
+            var holesfile=options.getfieldvalue('holes');
+            if ((typeof holesfile) === 'string') {
+                console.log('bamg error message: file ' + holesfile + ' loading from file path not supported - holes must be file object (.shp or .exp)');
+            } else if ((typeof holesfile) === 'object') {
+                holes = holesfile;
+            } else {
+                console.log('"holes" type not supported yet - ' + (typeof holesfile));
+            }
+        }
+        var subdomains = [];
+        if (options.exist('subdomains')) {
+            var subdomainsfile=options.getfieldvalue('subdomains');
+            if ((typeof subdomainsfile) === 'string') {
+                console.log('bamg error message: file ' + subdomainsfile + ' loading from file path not supported - subdomains must be file object (.shp or .exp)');
+            } else if ((typeof subdomainsfile) === 'object') {
+                subdomains = subdomainsfile;
+            } else {
+                console.log('"subdomains" type not supported yet - ' + (typeof subdomainsfile));
+            }
+        }
+
+        //Build geometry 
+        var count=0;
+        for (var i=0; i < domain.length; i++) {
+
+            //Check that the domain is closed
+            if (domain[i].x[0] != domain[i].x[domain[i].x.length-1] || domain[i].y[0] != domain[i].y[domain[i].y.length-1]) {
+                console.log('bamg error message: all contours provided in "domain" should be closed');
+            }
+
+            //TODO: Implement ContourToNodes
+            //Checks that all holes are INSIDE the principle domain outline
+            //if (i>1) {
+            //    flags=ContourToNodes(domain[i].x,domain[i].y,domain[0],0);
+            //    if (ArrayAny(ArrayFlip(flags))) {
+            //        console.log('bamg error message: All holes should be strictly inside the principal domain');
+            //    }
+            //}
+            console.log('bamg warning message: All holes should be strictly inside the principal domain. No checks are currently implemented.');
+
+            //Check orientation
+            var nods = domain[i].nods-1; //the domain are closed 1=end;
+            var test = ArraySum(ArrayAdd(ArraySubtract(domain[i].x.slice(1,nods), ArrayMultiply(domain[i].x.slice(0,nods-1), domain[i].y.slice(1,nods))), domain[i].y.slice(0,nods-1)));
+            if ((i==0 && test>0) || (i>0 && test<0)) {
+                console.log('At least one contour was not correctly oriented and has been re-oriented');
+                //TODO: Implement flipud
+                //domain[i].x = flipud(domain[i].x); domain[i].y = flipud(domain[i].y);
+            }
+
+            //Add all points to bamg_geometry
+            for (var j = 0; j < nods; j++) {
+                vertices.push(domain[i].x[j]);
+                vertices.push(domain[i].y[j]);
+                vertices.push(1);
+            }
+            for (var j = 0; j < nods; j++) {
+                edges.push(ArrayRange(count + 1, count + nods));
+                edges.push(ArrayConcat(ArrayRange(count + 2, count + nods), [count + 1]));
+                edges.push(1);
+            }
+            subdomains.push([2, count + 1, 1, -subdomain_ref]);
+            subdomain_ref = subdomain_ref + 1;
+            bamg_geometry.Vertices = vertices;
+            bamg_geometry.Edges = edges;
+            bamg_geometry.SubDomains = subdomains;
+
+            //update counter
+            count=count+nods;
+        }
+        for (var i=0; i < holes.length; i++) {
+
+            //Check that the subdomains is closed
+            if (holes[i].x[0] != holes[i].x[holes[i].x.length-1] || holes[i].y[0] != holes[i].y[holes[i].y.length-1]) {
+                console.log('bamg error message: all contours provided in "domain" should be closed');
+            }
+            //Checks that all holes are INSIDE the principle domain outline
+            //flags=ContourToNodes(holes[i].x,holes[i].y,domain[0],0);
+            //if ArrayAny(ArrayFlip(flags)) { console.log('bamg error message: All holes should be strictly inside the principal domain'); }
+
+            //TODO: Implement ContourToNodes
+            //Checks that all holes are INSIDE the principle domain outline
+            //if (i>1) {
+            //    flags=ContourToNodes(domain[i].x,domain[i].y,domain[0],0);
+            //    if (ArrayAny(ArrayFlip(flags))) {
+            //        console.log('bamg error message: All holes should be strictly inside the principal domain');
+            //    }
+            //}
+            console.log('bamg warning message: all holes should be strictly inside the principal domain. no checks are currently implemented.');
+
+            //Check that hole is correctly oriented
+            var nods = holes[i].nods-1; //the holes are closed 1=end;
+            var test = ArraySum(ArrayAdd(ArraySubtract(holes[i].x.slice(1,nods), ArrayMultiply(holes[i].x.slice(0,nods-1), holes[i].y.slice(1,nods))), holes[i].y.slice(0,nods-1)));
+            if ((i==0 && test>0) || (i>0 && test<0)) {
+                console.log('At least one contour was not correctly oriented and has been re-oriented');
+                //TODO: Implement flipud
+                //holes[i].x = flipud(holes[i].x); holes[i].y = flipud(holes[i].y);
+            }
+
+            //Add all points to bamg_geometry
+            for (var j = 0; j < nods; j++) {
+                vertices.push(holes[i].x[j]);
+                vertices.push(holes[i].y[j]);
+                vertices.push(1);
+            }
+            for (var j = 0; j < nods; j++) {
+                edges.push(ArrayRange(count + 1, count + nods));
+                edges.push(ArrayConcat(ArrayRange(count + 2, count + nods), [count + 1]));
+                edges.push(1);
+            }
+            subdomains.push([2, count + 1, 1, -hole_ref]);
+            hole_ref = hole_ref + 1;
+            bamg_geometry.Vertices = vertices;
+            bamg_geometry.Edges = edges;
+            bamg_geometry.SubDomains = subdomains;
+
+            //update counter
+            count=count+nods;
+        }
+        for (var i=0; i < subdomains.length; i++) {
+
+            //Check that the subdomains is closed
+            if (subdomains[i].x[0] != subdomains[i].x[subdomains[i].x.length-1] || subdomains[i].y[0] != subdomains[i].y[subdomains[i].y.length-1]) {
+                console.log('bamg error message: all contours provided in "subdomains" should be closed');
+            }
+
+            //TODO: Implement ContourToNodes
+            //Checks that all holes are INSIDE the principle domain outline
+            //flags=ContourToNodes(subdomains[i].x,subdomains[i].y,domain[0],0);
+            //if ArrayAny(ArrayFlip(flags)) {
+            //    console.log('bamg error message: All holes should be strictly inside the principal domain');
+            //}
+            console.log('bamg warning message: all holes should be strictly inside the principal domain. no checks are currently implemented.');
+
+            //Check that hole is correctly oriented
+            var nods=subdomains[i].nods-1; //the subdomains are closed 1=end;
+            var test = ArraySum(ArrayAdd(ArraySubtract(subdomains[i].x.slice(1,nods), ArrayMultiply(subdomains[i].x.slice(0,nods-1), subdomains[i].y.slice(1,nods))), subdomains[i].y.slice(0,nods-1)));
+            if ((i==0 && test>0) || (i>0 && test<0)) {
+                console.log('At least one contour was not correctly oriented and has been re-oriented');
+                //TODO: Implement flipud
+                //subdomains[i].x = flipud(subdomains[i].x); subdomains[i].y = flipud(subdomains[i].y);
+            }
+
+            //Add all points to bamg_geometry
+            for (var j = 0; j < nods; j++) {
+                vertices.push(subdomains[i].x[j]);
+                vertices.push(subdomains[i].y[j]);
+                vertices.push(1);
+            }
+            for (var j = 0; j < nods; j++) {
+                edges.push(ArrayRange(count + 1, count + nods));
+                edges.push(ArrayConcat(ArrayRange(count + 2, count + nods), [count + 1]));
+                edges.push(1);
+            }
+            subdomains.push([2, count + 1, 1, subdomain_ref]);
+            subdomain_ref = subdomain_ref + 1;
+            bamg_geometry.Vertices = vertices;
+            bamg_geometry.Edges = edges;
+            bamg_geometry.SubDomains = subdomains;
+
+            //update counter
+            count=count+nods;
+        }
+        if (options.getfieldvalue('vertical',0)) {
+            if (options.getfieldvalue('Markers',[]).length != bamg_geometry.Edges.length) {
+                console.log('for 2d vertical mesh, "Markers" option is required, and should be of size ' + bamg_geometry.Edges.length);
+            }
+        }
+        if (options.getfieldvalue('Markers',[]).length == bamg_geometry.Edges.length) {
+                console.log('BAMG Markers not implemented yet.');
+            //bamg_geometry.Edges(:,3)=options.getfieldvalue('Markers');
+        }
+        /*
+
+        //take care of rifts
+        if options.exist('rifts') {
+
+            //Check that file exists
+            riftfile=options.getfieldvalue('rifts');
+            [pathr,namer,extr]=fileparts(riftfile);
+            if !exist(riftfile {'file')
+                console.log(['bamg error message: file ' riftfile ' not found ']);
+            } else if strcmp(extr,'.exp') {
+                rift=expread(riftfile);
+            } else if strcmp(extr,'.shp') {
+                rift=shpread(riftfile);
+            }
+            //read rift file according to its extension: 
+            [path,name,ext]=fileparts(riftfile);
+            if strcmp(ext,'.exp') {
+                rift=expread(riftfile);
+            } else if strcmp(ext,'.shp') {
+                rift=shpread(riftfile);
+            } else {
+                console.log(['bamg error message: file ' riftfile ' format not supported (.shp or .exp)']);
+            }
+
+            for i=1:length(rift) {
+
+                //detect whether all points of the rift are inside the domain
+                flags=ContourToNodes(rift[i].x,rift[i].y,domain[0],0);
+                if (ArrayFlip(flags)) {
+                    console.log('one rift has all its points outside of the domain outline'),
+
+                } else if (ArrayAny(ArrayFlip(flags))) {
+                    //We LOTS of work to do
+                    console.log('Rift tip outside of or on the domain has been detected and is being processed...');
+
+                    //check that only one point is outside (for now)
+                    if (ArraySum(ArrayFlip(flags))!=1) {
+                        console.log('bamg error message: only one point outside of the domain is supported yet');
+                    }
+
+                    //Move tip outside to the first position
+                    if (flags[0]==0) {
+                        //OK, first point is outside (do nothing),
+                    } else if (flags[flags.length-1]==0) {
+                        rift[i].x=flipud(rift[i].x);
+                        rift[i].y=flipud(rift[i].y);
+                    } else {
+                        console.log('bamg error message: only a rift tip can be outside of the domain');
+                    }
+
+                    //Get cordinate of intersection point
+                    x1=rift[i].x[0]; y1=rift[i].y[0];
+                    x2=rift[i].x[1]; y2=rift[i].y[1];
+                    for (var j=0; j < domain[0].x)-1 j++) {
+                        if SegIntersect([x1 y1; x2 y2],[domain[0].x(j) domain[0].y(j); domain[0].x(j+1) domain[0].y(j+1)]) {
+
+                            //Get position of the two nodes of the edge in domain
+                            i1=j;
+                            i2=j+1;
+
+                            //rift is crossing edge [i1 i2] of the domain
+                            //Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
+                            x3=domain[0].x[i1]; y3=domain[0].y[i1];
+                            x4=domain[0].x[i2]; y4=domain[0].y[i2];
+                            x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+                            y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+
+                            segdis= sqrt((x4-x3)^2+(y4-y3)^2);
+                            tipdis=[sqrt((x-x3)^2+(y-y3)^2)  sqrt((x-x4)^2+(y-y4)^2)];
+
+                            if (min(tipdis)/segdis) < options.getfieldvalue('toltip',0) {
+                                disp('moving tip-domain intersection point');
+
+                                //Get position of the closer point
+                                if tipdis[0]>tipdis[1] {
+                                    pos=i2;
+                                } else {
+                                    pos=i1;
+                                }
+
+                                //This point is only in Vertices (number pos).
+                                //OK, now we can add our own rift
+                                nods=rift[i].nods-1;
+                                bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift[i].x(2:end) rift[i].y(2:end) ones(nods,1)]];
+                                bamg_geometry.Edges=[bamg_geometry.Edges;...
+                                    pos count+1  (1+i);...
+                                    [transpose(count+1:count+nods-1) transpose(count+2:count+nods)  (1+i)*ones(nods-1,1)]];
+                                count=count+nods;
+
+                                break;
+
+                            } else {
+                                //Add intersection point to Vertices
+                                bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1];
+                                count=count+1;
+
+                                //Decompose the crossing edge into 2 subedges
+                                pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2);
+                                if isempty(pos) console.log('bamg error message: a problem occurred...'); }
+                                bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);...
+                                    bamg_geometry.Edges(pos,1) count                      bamg_geometry.Edges(pos,3);...
+                                    count                      bamg_geometry.Edges(pos,2) bamg_geometry.Edges(pos,3);...
+                                    bamg_geometry.Edges(pos+1:end,:)];
+
+                                //OK, now we can add our own rift
+                                nods=rift[i].nods-1;
+                                bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift[i].x(2:end) rift[i].y(2:end) ones(nods,1)]];
+                                bamg_geometry.Edges=[bamg_geometry.Edges;...
+                                    count  count+1  2 ;...
+                                    [transpose(count+1:count+nods-1) transpose(count+2:count+nods)  (1+i)*ones(nods-1,1)]];
+                                count=count+nods;
+
+                                break;
+                            }
+                        }
+                    }
+                } else {
+                    nods=rift[i].nods-1;
+                    bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift[i].x(:) rift[i].y(:) ones(nods+1,1)]];
+                    bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose(count+2:count+nods+1)  (1+i)*ones(nods,1)]];
+                    count=count+nods+1;
+                }
+            }
+        }
+
+        //Deal with tracks
+        if options.exist('tracks') {
+
+            //read tracks
+            track=options.getfieldvalue('tracks');
+            if all(ischar(track)) {
+                A=expread(track);
+                track=[];
+                for i=1:length(A), 
+                    track=[track; [A[i].x A[i].y]];
+                }
+            } else {
+                track=double(track); //for some reason, it is of class "single"
+            }
+            if(size(track,2)==2), track=[track 3.*ones(size(track,1),1)]; }
+
+            //only keep those inside
+            flags=ContourToNodes(track(:,1),track(:,2),domainfile,0);
+            track=track(find(flags),:);
+
+            //Add all points to bamg_geometry
+            nods=size(track,1);
+            bamg_geometry.Vertices=[bamg_geometry.Vertices; track];
+            bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose(count+2:count+nods)  3.*ones(nods-1,1)]];
+
+            //update counter
+            count=count+nods;
+        }
+
+        //Deal with vertices that need to be kept by mesher
+        if (options.exist('RequiredVertices')) {
+
+            //recover RequiredVertices
+            requiredvertices= options.getfieldvalue('RequiredVertices'); //for some reason, it is of class "single"
+            if (requiredvertices[0].length==2) {
+                 requiredvertices=[requiredvertices 4.*ones(size(requiredvertices,1),1)];
+            }    
+
+            //only keep those inside
+            flags=ContourToNodes(requiredvertices(:,1),requiredvertices(:,2),domain[0],0);
+            requiredvertices=requiredvertices(find(flags),:);
+
+            //Add all points to bamg_geometry
+            nods=size(requiredvertices,1);
+            bamg_geometry.Vertices=[bamg_geometry.Vertices; requiredvertices];
+
+            //update counter
+            count=count+nods;
+
+        }
+
+        //Deal with RequiredEdges
+        if (options.getfieldvalue('NoBoundaryRefinment',0)==1) {
+            bamg_geometry.RequiredEdges=[1:size(bamg_geometry.Edges,1)]';
+        }
+
+        //process geom
+        //bamg_geometry=processgeometry(bamg_geometry,options.getfieldvalue('tol',NaN),domain[0]);
+        */
+
+    } /*else if ((typeof md.private.bamg === 'object') && ('geometry' in md.private.bamg)) {
+        bamg_geometry=bamggeom(md.private.bamg.geometry); 
+    } else {
+        //do nothing...
+    }
+    //}}}
+    // Bamg Mesh parameters {{{
+    if (!options.exist('domain') && md.mesh.numberofvertices!=0 && md.mesh.elementtype()=='Tria') {
+
+        if ((typeof md.private.bamg === 'object') && ('mesh' in md.private.bamg)) {
+            bamg_mesh=bamgmesh(md.private.bamg.mesh);
+        } else {
+            bamg_mesh.Vertices=[md.mesh.x md.mesh.y ones(md.mesh.numberofvertices,1)];
+            bamg_mesh.Triangles=[md.mesh.elements ones(md.mesh.numberofelements,1)];
+        }
+
+        if (typeof md.rifts.riftstruct === 'object') {
+            console.log('bamg error message: rifts not supported yet. Do meshprocessrift AFTER bamg');
+        }
+    }*/
+    //}}}
+    // Bamg Options {{{
+    bamg_options.Crack=options.getfieldvalue('Crack',0);
+    bamg_options.anisomax=options.getfieldvalue('anisomax',Math.pow(10,30));
+    bamg_options.coeff=options.getfieldvalue('coeff',1.);
+    bamg_options.cutoff=options.getfieldvalue('cutoff',Math.pow(10,-5));
+    bamg_options.err=options.getfieldvalue('err',0.01);
+    bamg_options.errg=options.getfieldvalue('errg',0.1);
+    bamg_options.field=options.getfieldvalue('field',[]);
+    bamg_options.gradation=options.getfieldvalue('gradation',1.5);
+    bamg_options.Hessiantype=options.getfieldvalue('Hessiantype',0);
+    bamg_options.hmin=options.getfieldvalue('hmin',Math.pow(10,-100));
+    bamg_options.hmax=options.getfieldvalue('hmax',Math.pow(10,100));
+    bamg_options.hminVertices=options.getfieldvalue('hminVertices',[]);
+    bamg_options.hmaxVertices=options.getfieldvalue('hmaxVertices',[]);
+    bamg_options.hVertices=options.getfieldvalue('hVertices',[]);
+    bamg_options.KeepVertices=options.getfieldvalue('KeepVertices',1);
+    bamg_options.maxnbv=options.getfieldvalue('maxnbv',Math.pow(10,6));
+    bamg_options.maxsubdiv=options.getfieldvalue('maxsubdiv',10.);
+    bamg_options.metric=options.getfieldvalue('metric',[]);
+    bamg_options.Metrictype=options.getfieldvalue('Metrictype',0);
+    bamg_options.nbjacobi=options.getfieldvalue('nbjacobi',1);
+    bamg_options.nbsmooth=options.getfieldvalue('nbsmooth',3);
+    bamg_options.omega=options.getfieldvalue('omega',1.8);
+    bamg_options.power=options.getfieldvalue('power',1.);
+    bamg_options.splitcorners=options.getfieldvalue('splitcorners',1);
+    bamg_options.verbose=options.getfieldvalue('verbose',1);
+    //}}}
+
+    //call Bamg
+    console.log("calling BamgMesher");
+    var return_array=BamgMesher(bamg_mesh,bamg_geometry,bamg_options);
+    var bamgmesh_out=return_array[0];
+    var bamggeom_out=return_array[1];
+//
+//    if (options.getfieldvalue('vertical',0)) {
+//        md.mesh=mesh2dvertical();
+//        md.mesh.x=bamgmesh_out.Vertices(:,1);
+//        md.mesh.y=bamgmesh_out.Vertices(:,2);
+//        md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
+//        md.mesh.edges=bamgmesh_out.IssmEdges;
+//        md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3);
+//        md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4);
+//
+//        //Fill in rest of fields:
+//        md.mesh.numberofelements=size(md.mesh.elements,1);
+//        md.mesh.numberofvertices=length(md.mesh.x);
+//        md.mesh.numberofedges=size(md.mesh.edges,1);
+//        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+//
+//        //Now, build the connectivity tables for this mesh.
+//        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+//
+//    } else if options.getfieldvalue('3dsurface',0) {
+//        
+//        md.mesh=mesh3dsurface();
+//        md.mesh.x=bamgmesh_out.Vertices(:,1);
+//        md.mesh.y=bamgmesh_out.Vertices(:,2);
+//        md.mesh.z=md.mesh.x; md.mesh.z(:)=0;
+//        md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
+//        md.mesh.edges=bamgmesh_out.IssmEdges;
+//        md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3);
+//        md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4);
+//
+//        //Fill in rest of fields:
+//        md.mesh.numberofelements=size(md.mesh.elements,1);
+//        md.mesh.numberofvertices=length(md.mesh.x);
+//        md.mesh.numberofedges=size(md.mesh.edges,1);
+//        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+//
+//    } else { 
+//        md.mesh=mesh2d();
+//        md.mesh.x=bamgmesh_out.Vertices(:,1);
+//        md.mesh.y=bamgmesh_out.Vertices(:,2);
+//        md.mesh.elements=bamgmesh_out.Triangles(:,1:3);
+//        md.mesh.edges=bamgmesh_out.IssmEdges;
+//        md.mesh.segments=bamgmesh_out.IssmSegments(:,1:3);
+//        md.mesh.segmentmarkers=bamgmesh_out.IssmSegments(:,4);
+//
+//        //Fill in rest of fields:
+//        md.mesh.numberofelements=size(md.mesh.elements,1);
+//        md.mesh.numberofvertices=length(md.mesh.x);
+//        md.mesh.numberofedges=size(md.mesh.edges,1);
+//        md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
+//    }
+//
+//    //Bamg private fields
+//    md.private.bamg=struct();
+//    md.private.bamg.mesh=bamgmesh(bamgmesh_out);
+//    md.private.bamg.geometry=bamggeom(bamggeom_out);
+//    md.mesh.elementconnectivity=md.private.bamg.mesh.ElementConnectivity;
+//    md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;
+//
+//    //Check for orphan
+//    if ArrayAny(!ismember(1:md.mesh.numberofvertices,sort(unique(reshape(md.mesh.elements,3*md.mesh.numberofelements,1)))))
+//        console.log('Output mesh has orphans. Check your Domain and/or RequiredVertices');
+//    }
+//    } 
+//
+//    function geom=processgeometry(geom,tol,outline) // {{{
+//
+//    //Deal with edges
+//    disp('Checking Edge crossing...');
+//    i=0;
+//    while (i<size(geom.Edges,1)),
+//
+//        //edge counter
+//        i=i+1;
+//
+//        //Get coordinates
+//        x1=geom.Vertices(geom.Edges(i,1),1);
+//        y1=geom.Vertices(geom.Edges(i,1),2);
+//        x2=geom.Vertices(geom.Edges(i,2),1);
+//        y2=geom.Vertices(geom.Edges(i,2),2);
+//        color1=geom.Edges(i,3);
+//
+//        j=i; //test edges located AFTER i only
+//        while (j<size(geom.Edges,1)),
+//
+//            //edge counter
+//            j=j+1;
+//
+//            //Skip if the two edges already have a vertex in common
+//            if ArrayAny(ismember(geom.Edges(i,1:2),geom.Edges(j,1:2))),
+//                continue
+//            }
+//
+//            //Get coordinates
+//            x3=geom.Vertices(geom.Edges(j,1),1);
+//            y3=geom.Vertices(geom.Edges(j,1),2);
+//            x4=geom.Vertices(geom.Edges(j,2),1);
+//            y4=geom.Vertices(geom.Edges(j,2),2);
+//            color2=geom.Edges(j,3);
+//
+//            //Check if the two edges are crossing one another
+//            if SegIntersect([x1 y1; x2 y2],[x3 y3; x4 y4]),
+//
+//                //Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html)
+//                x=det([det([x1 y1; x2 y2])  x1-x2;det([x3 y3; x4 y4])  x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+//                y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
+//
+//                //Add vertex to the list of vertices
+//                geom.Vertices(end+1,:)=[x y min(color1,color2)];
+//                id=size(geom.Vertices,1);
+//
+//                //Update edges i and j
+//                edgei=geom.Edges(i,:);
+//                edgej=geom.Edges(j,:);
+//                geom.Edges(i,:)    =[edgei[0] id       edgei(3)];
+//                geom.Edges(end+1,:)=[id       edgei[1] edgei(3)];
+//                geom.Edges(j,:)    =[edgej[0] id       edgej(3)];
+//                geom.Edges(end+1,:)=[id       edgej[1] edgej(3)];
+//
+//                //update current edge second tip
+//                x2=x; y2=y;
+//            }
+//        }
+//
+//    }
+//
+//    //Check point outside
+//    disp('Checking for points outside the domain...');
+//    i=0;
+//    num=0;
+//    while (i<size(geom.Vertices,1)),
+//
+//        //vertex counter
+//        i=i+1;
+//
+//        //Get coordinates
+//        x=geom.Vertices(i,1);
+//        y=geom.Vertices(i,2);
+//        color=geom.Vertices(i,3);
+//
+//        //Check that the point is inside the domain
+//        if (color!=1 & !ContourToNodes(x,y,outline[0],1)),
+//
+//            //Remove points from list of Vertices
+//            num=num+1;
+//            geom.Vertices(i,:)=[];
+//
+//            //update edges
+//            [posedges dummy]=find(geom.Edges==i);
+//            geom.Edges(posedges,:)=[];
+//            posedges=find(geom.Edges>i);
+//            geom.Edges(posedges)=geom.Edges(posedges)-1;
+//
+//            //update counter
+//            i=i-1;
+//        }
+//    }
+//    if num,
+//        disp(['WARNING: ' num2str(num) ' points outside the domain outline have been removed']);
+//    }
+//
+//    //Check point spacing
+//    if !isnan(tol),
+//        disp('Checking point spacing...');
+//        i=0;
+//        while (i<size(geom.Vertices,1)),
+//
+//            //vertex counter
+//            i=i+1;
+//
+//            //Get coordinates
+//            x1=geom.Vertices(i,1);
+//            y1=geom.Vertices(i,2);
+//
+//            j=i; //test edges located AFTER i only
+//            while (j<size(geom.Vertices,1)),
+//
+//                //vertex counter
+//                j=j+1;
+//
+//                //Get coordinates
+//                x2=geom.Vertices(j,1);
+//                y2=geom.Vertices(j,2);
+//
+//                //Check whether the two vertices are too close
+//                if ((x2-x1)^2+(y2-y1)^2<tol^2)
+//
+//                    //Remove points from list of Vertices
+//                    geom.Vertices(j,:)=[];
+//
+//                    //update edges
+//                    posedges=find(ismember(geom.Edges,j));
+//                    geom.Edges(posedges)=i;
+//                    posedges=find(geom.Edges>j);
+//                    geom.Edges(posedges)=geom.Edges(posedges)-1;
+//
+//                    //update counter
+//                    j=j-1;
+//
+//                }
+//            }
+//        }
+//    }
+//    //remove empty edges
+//    geom.Edges(find(geom.Edges(:,1)==geom.Edges(:,2)),:)=[];
+} // }}}
Index: /issm/trunk-jpl/src/wrappers/BamgMesher/BamgMesher.h
===================================================================
--- /issm/trunk-jpl/src/wrappers/BamgMesher/BamgMesher.h	(revision 22876)
+++ /issm/trunk-jpl/src/wrappers/BamgMesher/BamgMesher.h	(revision 22877)
@@ -45,4 +45,16 @@
 #endif
 
+#ifdef _HAVE_JAVASCRIPT_MODULES_
+/* serial input macros: */
+#define BAMGMESHIN Vertices_mesh_in,VerticesSize_mesh_in,Edges_mesh_in,EdgesSize_mesh_in,Triangles_mesh_in,TrianglesSize_mesh_in,CrackedEdges_mesh_in,CrackedEdgesSize_mesh_in,VerticesOnGeomEdge_mesh_in,VerticesOnGeomEdgeSize_mesh_in,VerticesOnGeomVertex_mesh_in,VerticesOnGeomVertexSize_mesh_in,EdgesOnGeomEdge_mesh_in,EdgesOnGeomEdgeSize_mesh_in,IssmSegments_mesh_in,IssmSegmentsSize_mesh_in
+#define BAMGGEOMIN Vertices_geom_in,VerticesSize_geom_in,Edges_geom_in,EdgesSize_geom_in,Corners_geom_in,CornersSize_geom_in,RequiredVertices_geom_in,RequiredVerticesSize_geom_in,RequiredCrackedEdges_geom_in,RequiredCrackedEdgesSize_geom_in,CrackedEdges_geom_in,CrackedEdgesSize_geom_in,SubDomains_geom_in,SubDomainsSize_geom_in
+#define BAMGOPTIONS pbamgopts,anisomax,cutoff,coeff,errg,gradation,Hessiantype,maxnbv,maxsubdiv,Metrictype,nbjacobi,nbsmooth,omega,power,verbose,Crack,KeepVertices,splitcorners,hmin,hmax,hminVertices,hmaxVertices,hVertices,metric,field,err
+/* serial output macros: */
+#define BAMGMESHOUT Vertices_mesh_out,VerticesSize_mesh_out,Edges_mesh_out,EdgesSize_mesh_out,Triangles_mesh_out,TrianglesSize_mesh_out,CrackedEdges_mesh_out,CrackedEdgesSize_mesh_out,VerticesOnGeomEdge_mesh_out,VerticesOnGeomEdgeSize_mesh_out,VerticesOnGeomVertex_mesh_out,VerticesOnGeomVertexSize_mesh_out,EdgesOnGeomEdge_mesh_out,EdgesOnGeomEdgeSize_mesh_out,IssmSegments_mesh_out,IssmSegmentsSize_mesh_out
+#define BAMGGEOMOUT Vertices_geom_out,VerticesSize_geom_out,Edges_geom_out,EdgesSize_geom_out,Corners_geom_out,CornersSize_geom_out,RequiredVertices_geom_out,RequiredVerticesSize_geom_out,RequiredCrackedEdges_geom_out,RequiredCrackedEdgesSize_geom_out,CrackedEdges_geom_out,CrackedEdgesSize_geom_out,SubDomains_geom_out,SubDomainsSize_geom_out
+#define WRAPPER(modulename) extern "C" { int  BamgMesherModule(double* Vertices_mesh_out, int* VerticesSize_mesh_out, double* Edges_mesh_out, int* EdgesSize_mesh_out, double* Triangles_mesh_out, int* TrianglesSize_mesh_out, double* CrackedEdges_mesh_out, int* CrackedEdgesSize_mesh_out, double* VerticesOnGeomEdge_mesh_out, int* VerticesOnGeomEdgeSize_mesh_out, double* VerticesOnGeomVertex_mesh_out, int* VerticesOnGeomVertexSize_mesh_out, double* EdgesOnGeomEdge_mesh_out, int* EdgesOnGeomEdgeSize_mesh_out, double* IssmSegments_mesh_out, int* IssmSegmentsSize_mesh_out, double* Vertices_geom_out, int* VerticesSize_geom_out, double* Edges_geom_out, int* EdgesSize_geom_out, double* Corners_geom_out, int* CornersSize_geom_out, double* RequiredVertices_geom_out, int* RequiredVerticesSize_geom_out, double* RequiredCrackedEdges_geom_out, int* RequiredCrackedEdgesSize_geom_out, double* CrackedEdges_geom_out, int* CrackedEdgesSize_geom_out, double* SubDomains_geom_out, int* SubDomainsSize_geom_out, double* Vertices_mesh_in, int* VerticesSize_mesh_in, double* Edges_mesh_in, int* EdgesSize_mesh_in, double* Triangles_mesh_in, int* TrianglesSize_mesh_in, double* CrackedEdges_mesh_in, int* CrackedEdgesSize_mesh_in, double* VerticesOnGeomEdge_mesh_in, int* VerticesOnGeomEdgeSize_mesh_in, double* VerticesOnGeomVertex_mesh_in, int* VerticesOnGeomVertexSize_mesh_in, double* EdgesOnGeomEdge_mesh_in, int* EdgesOnGeomEdgeSize_mesh_in, double* IssmSegments_mesh_in, int* IssmSegmentsSize_mesh_in, double* Vertices_geom_in, int* VerticesSize_geom_in, double* Edges_geom_in, int* EdgesSize_geom_in, double* Corners_geom_in, int* CornersSize_geom_in, double* RequiredVertices_geom_in, int* RequiredVerticesSize_geom_in, double* RequiredCrackedEdges_geom_in, int* RequiredCrackedEdgesSize_geom_in, double* CrackedEdges_geom_in, int* CrackedEdgesSize_geom_in, double* SubDomains_geom_in, int* SubDomainsSize_geom_in, double anisomax, double cutoff, double coeff, double errg, double gradation, int Hessiantype, int maxnbv, double maxsubdiv, int Metrictype, int nbjacobi, int nbsmooth, double omega, double power, int verbose, int Crack, int KeepVertices, int splitcorners, double hmin, double hmax, double* hminVertices, double* hmaxVertices, double* hVertices, double* metric, double* field, double err)
+#define _DO_NOT_LOAD_GLOBALS_//we only load globals for TriangleModule.js, not other modules!
+#endif
+
 /* serial arg counts: */
 #undef NLHS
Index: /issm/trunk-jpl/src/wrappers/BamgMesher/BamgMesher.js
===================================================================
--- /issm/trunk-jpl/src/wrappers/BamgMesher/BamgMesher.js	(revision 22877)
+++ /issm/trunk-jpl/src/wrappers/BamgMesher/BamgMesher.js	(revision 22877)
@@ -0,0 +1,272 @@
+function doublesToHeap(array) {
+	var doubleArray = new Float64Array(array);
+    var numBytes = doubleArray.length * doubleArray.BYTES_PER_ELEMENT;
+	var doubleArrayPtr = Module._malloc(numBytes);
+    var doubleArrayHeap = new Uint8Array(Module.HEAPU8.buffer, doubleArrayPtr, numBytes);
+	doubleArrayHeap.set(new Uint8Array(doubleArray.buffer));
+    return doubleArrayHeap.byteOffset;
+}
+function intsToHeap(array) {
+	var intArray = new Int32Array(array);
+    var numBytes = intArray.length * intArray.BYTES_PER_ELEMENT;
+	var intArrayPtr = Module._malloc(numBytes);
+    var intArrayHeap = new Uint8Array(Module.HEAPU8.buffer, intArrayPtr, numBytes);
+	intArrayHeap.set(new Uint8Array(intArray.buffer));
+    return intArrayHeap.byteOffset;
+}
+function heapToDoubles(pptr, size) {
+	var ptr = Module.getValue(pptr,'i32');
+	var array = Module.HEAPF64.slice(ptr / 8, ptr / 8 + size[0] * size[1]);
+	return ListToMatrix(array, size[1]);
+}
+function heapToInts(pptr, nods) {
+	var ptr = Module.getValue(pptr,'i32');
+	return Module.HEAPU32.slice(ptr / 4, ptr / 4 + nods);
+}
+function BamgMesher(md, bamgmesh, bamggeom, bamgopts) {
+/*
+	   usage: var array = Triangle(domain,rifts,area);
+	      where: array is made of [index,x,y,segments,segmentmarkers]
+		  and index,x,y defines a triangulation, segments is an array made 
+	      of exterior segments to the mesh domain outline, segmentmarkers is an array 
+		  flagging each segment, domain a js array defining the domain outline  (sames for 
+		  rifts) and area is the maximum area desired for any element of the resulting mesh.
+
+		  Ok, for now, we are not dealing with rifts. Also, the domain is made of only one 
+		  profile, this to avoid passing a double** pointer to js. 
+*/
+
+	//Dynamic allocations: {{{
+	//Retrieve domain arrays, and allocate on Module heap: 
+	//input
+    var pVertices_mesh_in                   = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.Vertices));
+    var pVerticesSize_mesh_in               = intsToHeap([bamgmesh.Vertices.length, bamgmesh.Vertices[0].length]);
+    var pEdges_mesh_in                      = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.Edges));
+    var pEdgesSize_mesh_in                  = intsToHeap([bamgmesh.Edges.length, bamgmesh.Edges[0].length]);
+    var pTriangles_mesh_in                  = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.Triangle));
+    var pTrianglesSize_mesh_in              = intsToHeap([bamgmesh.Triangle.length, bamgmesh.Triangle[0].length]);
+    var pCrackedEdges_mesh_in               = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.CrackedEdges));
+    var pCrackedEdgesSize_mesh_in           = intsToHeap([bamgmesh.CrackedEdges.length, bamgmesh.CrackedEdges[0].length]);
+    var pVerticesOnGeomEdge_mesh_in         = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.VerticesOnGeomEdge));
+    var pVerticesOnGeomEdgeSize_mesh_in     = intsToHeap([bamgmesh.VerticesOnGeomEdge.length, bamgmesh.VerticesOnGeomEdge[0].length]);
+    var pVerticesOnGeomVertex_mesh_in       = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.VerticesOnGeomVertex));
+    var pVerticesOnGeomVertexSize_mesh_in   = intsToHeap([bamgmesh.VerticesOnGeomVertex.length, bamgmesh.VerticesOnGeomVertex[0].length]);
+    var pEdgesOnGeomEdge_mesh_in            = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.EdgesOnGeomEdge));
+    var pEdgesOnGeomEdgeSize_mesh_in        = intsToHeap([bamgmesh.EdgesOnGeomEdge.length, bamgmesh.EdgesOnGeomEdge[0].length]);
+    var pIssmSegments_mesh_in               = doublesToHeap(Array.prototype.concat.apply([], bamgmesh.IssmSegments));
+    var pIssmSegmentsSize_mesh_in           = intsToHeap([bamgmesh.IssmSegments.length, bamgmesh.IssmSegments[0].length]);
+
+    var pVertices_geom_in                   = doublesToHeap(Array.prototype.concat.apply([], bamggeom.Vertices));
+    var pVerticesSize_geom_in               = intsToHeap([bamggeom.Vertices.length, bamggeom.Vertices[0].length]);
+    var pEdges_geom_in                      = doublesToHeap(Array.prototype.concat.apply([], bamggeom.Edges));
+    var pEdgesSize_geom_in                  = intsToHeap([bamggeom.Edges.length, bamggeom.Edges[0].length]);
+    var pCorners_geom_in                    = doublesToHeap(Array.prototype.concat.apply([], bamggeom.Corners));
+    var pCornersSize_geom_in                = intsToHeap([bamggeom.Corners.length, bamggeom.Corners[0].length]);
+    var pRequiredVertices_geom_in           = doublesToHeap(Array.prototype.concat.apply([], bamggeom.RequiredVertices));
+    var pRequiredVerticesSize_geom_in       = intsToHeap([bamggeom.RequiredVertices.length, bamggeom.RequiredVertices[0].length]);
+    var pRequiredCrackedEdges_geom_in       = doublesToHeap(Array.prototype.concat.apply([], bamggeom.RequiredCrackedEdges));
+    var pRequiredCrackedEdgesSize_geom_in   = intsToHeap([bamggeom.RequiredCrackedEdges.length, bamggeom.RequiredCrackedEdges[0].length]);
+    var pCrackedEdges_geom_in               = doublesToHeap(Array.prototype.concat.apply([], bamggeom.CrackedEdges));
+    var pCrackedEdgesSize_geom_in           = intsToHeap([bamggeom.CrackedEdges.length, bamggeom.CrackedEdges[0].length]);
+    var pSubDomains_geom_in                 = doublesToHeap(Array.prototype.concat.apply([], bamggeom.SubDomains));
+    var pSubDomainsSize_geom_in             = intsToHeap([bamggeom.SubDomains.length, bamggeom.SubDomains[0].length]);
+
+    var anisomax                            = bamgopts.anisomax;
+    var cutoff                              = bamgopts.coeff;
+    var coeff                               = bamgopts.cutoff;
+    var errg                                = bamgopts.errg;
+    var gradation                           = bamgopts.gradation;
+    var Hessiantype                         = bamgopts.Hessiantype;
+    var maxnbv                              = bamgopts.maxnbv;
+    var maxsubdiv                           = bamgopts.maxsubdiv;
+    var Metrictype                          = bamgopts.Metrictype;
+    var nbjacobi                            = bamgopts.nbjacobi;
+    var nbsmooth                            = bamgopts.nbsmooth;
+    var omega                               = bamgopts.omega;
+    var power                               = bamgopts.power;
+    var verbose                             = bamgopts.verbose;
+    var Crack                               = bamgopts.Crack;
+    var KeepVertices                        = bamgopts.KeepVertices;
+    var splitcorners                        = bamgopts.splitcorners;
+    var hmin                                = bamgopts.hmin;
+    var hmax                                = bamgopts.hmax;
+    var phminVertices                       = doublesToHeap(Array.prototype.concat.apply([], bamgopts.hminVertices));
+    var phminVerticesSize                   = intsToHeap([bamgopts.hminVertices.length, bamgopts.hminVertices[0].length]);
+    var phmaxVertices                       = doublesToHeap(Array.prototype.concat.apply([], bamgopts.hmaxVertices));
+    var phmaxVerticesSize                   = intsToHeap([bamgopts.hmaxVertices.length, bamgopts.hmaxVertices[0].length]);
+    var phVertices                          = doublesToHeap(Array.prototype.concat.apply([], bamgopts.hVertices));
+    var phVerticesSize                      = intsToHeap([bamgopts.hVertices.length, bamgopts.hVertices[0].length]);
+    var pmetric                             = doublesToHeap(Array.prototype.concat.apply([], bamgopts.metric));
+    var pmetricSize                         = intsToHeap([bamgopts.metric.length, bamgopts.metric[0].length]);
+    var pfield                              = doublesToHeap(Array.prototype.concat.apply([], bamgopts.field));
+    var pfieldSize                          = intsToHeap([bamgopts.field.length, bamgopts.field[0].length]);
+    var err                                 = bamgopts.err;
+	
+	//output
+    var pVertices_geom_out                  = Module._malloc(4);
+    var pVerticesSize_geom_out              = Module._malloc(4); 
+    var pEdges_geom_out                     = Module._malloc(4); 
+    var pEdgesSize_geom_out                 = Module._malloc(4); 
+    var pCorners_geom_out                   = Module._malloc(4); 
+    var pCornersSize_geom_out               = Module._malloc(4); 
+    var pRequiredVertices_geom_out          = Module._malloc(4); 
+    var pRequiredVerticesSize_geom_out      = Module._malloc(4); 
+    var pRequiredCrackedEdges_geom_out      = Module._malloc(4); 
+    var pRequiredCrackedEdgesSize_geom_out  = Module._malloc(4); 
+    var pCrackedEdges_geom_out              = Module._malloc(4); 
+    var pCrackedEdgesSize_geom_out          = Module._malloc(4); 
+    var pSubDomains_geom_out                = Module._malloc(4); 
+    var pSubDomainsSize_geom_out            = Module._malloc(4); 
+
+    var pVertices_mesh_out                  = Module._malloc(4); 
+    var pVerticesSize_mesh_out              = Module._malloc(4); 
+    var pEdges_mesh_out                     = Module._malloc(4); 
+    var pEdgesSize_mesh_out                 = Module._malloc(4); 
+    var pTriangles_mesh_out                 = Module._malloc(4); 
+    var pTrianglesSize_mesh_out             = Module._malloc(4); 
+    var pCrackedEdges_mesh_out              = Module._malloc(4); 
+    var pCrackedEdgesSize_mesh_out          = Module._malloc(4); 
+    var pVerticesOnGeomEdge_mesh_out        = Module._malloc(4); 
+    var pVerticesOnGeomEdgeSize_mesh_out    = Module._malloc(4); 
+    var pVerticesOnGeomVertex_mesh_out      = Module._malloc(4); 
+    var pVerticesOnGeomVertexSize_mesh_out  = Module._malloc(4); 
+    var pEdgesOnGeomEdge_mesh_out           = Module._malloc(4); 
+    var pEdgesOnGeomEdgeSize_mesh_out       = Module._malloc(4); 
+    var pIssmSegments_mesh_out              = Module._malloc(4); 
+    var pIssmSegmentsSize_mesh_out          = Module._malloc(4); 
+	//}}}
+
+	//Declare BamgMesher module: 
+	BamgMesherModule = Module.cwrap('BamgMesherModule', 'number',[
+        'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 
+        'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 
+        'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 
+        'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 
+        'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number', 'number']);
+	
+	//Call BamgMesher module: 
+	BamgMesherModule(
+        pVertices_mesh_out, pVerticesSize_mesh_out, pEdges_mesh_out, pEdgesSize_mesh_out, pTriangles_mesh_out, pTrianglesSize_mesh_out, pCrackedEdges_mesh_out, pCrackedEdgesSize_mesh_out, pVerticesOnGeomEdge_mesh_out, pVerticesOnGeomEdgeSize_mesh_out, pVerticesOnGeomVertex_mesh_out, pVerticesOnGeomVertexSize_mesh_out, pEdgesOnGeomEdge_mesh_out, pEdgesOnGeomEdgeSize_mesh_out, pIssmSegments_mesh_out, pIssmSegmentsSize_mesh_out, 
+        pVertices_geom_out, pVerticesSize_geom_out, pEdges_geom_out, pEdgesSize_geom_out, pCorners_geom_out, pCornersSize_geom_out, pRequiredVertices_geom_out, pRequiredVerticesSize_geom_out, pRequiredCrackedEdges_geom_out, pRequiredCrackedEdgesSize_geom_out, pCrackedEdges_geom_out, pCrackedEdgesSize_geom_out, pSubDomaouts_geom_out, pSubDomaoutsSize_geom_out, 
+        pVertices_mesh_in, pVerticesSize_mesh_in, pEdges_mesh_in, pEdgesSize_mesh_in, pTriangles_mesh_in, pTrianglesSize_mesh_in, pCrackedEdges_mesh_in, pCrackedEdgesSize_mesh_in, pVerticesOnGeomEdge_mesh_in, pVerticesOnGeomEdgeSize_mesh_in, pVerticesOnGeomVertex_mesh_in, pVerticesOnGeomVertexSize_mesh_in, pEdgesOnGeomEdge_mesh_in, pEdgesOnGeomEdgeSize_mesh_in, pIssmSegments_mesh_in, pIssmSegmentsSize_mesh_in, 
+        pVertices_geom_in, pVerticesSize_geom_in, pEdges_geom_in, pEdgesSize_geom_in, pCorners_geom_in, pCornersSize_geom_in, pRequiredVertices_geom_in, pRequiredVerticesSize_geom_in, pRequiredCrackedEdges_geom_in, pRequiredCrackedEdgesSize_geom_in, pCrackedEdges_geom_in, pCrackedEdgesSize_geom_in, pSubDomains_geom_in, pSubDomainsSize_geom_in, 
+        anisomax, cutoff, coeff, errg, gradation, Hessiantype, maxnbv, maxsubdiv, Metrictype, nbjacobi, nbsmooth, omega, power, verbose, Crack, KeepVertices, splitcorners, hmin, hmax, phminVertices, phminVerticesSize, phmaxVertices, phmaxVerticesSize, phVertices, phVerticesSize, pmetric, pmetricSize, pfield, pfieldSize, err);
+	
+	/*Dynamic copying from heap: {{{*/
+	//recover mesh: 
+    var bamgmeshout = bamgmesh();
+    var bamggeomout = bamggeom();
+
+    bamgmeshout.VerticesSize                = heapToint(pVerticesSize_mesh_out, 2);
+    bamgmeshout.Vertices                    = heapToDoubles(heapToDoubles(pVertices_mesh_out), bamgmeshout.VerticesSize);
+    bamgmeshout.EdgesSize                   = heapToint(pEdgesSize_mesh_out, 2);
+    bamgmeshout.Edges                       = heapToDoubles(heapToDoubles(pEdges_mesh_out), bamgmeshout.EdgesSize);
+    bamgmeshout.TrianglesSize               = heapToint(pTrianglesSize_mesh_out, 2);
+    bamgmeshout.Triangles                   = heapToDoubles(heapToDoubles(pTriangles_mesh_out), bamgmeshout.TrianglesSize);
+    bamgmeshout.CrackedEdgesSize            = heapToint(pCrackedEdgesSize_mesh_out, 2);
+    bamgmeshout.CrackedEdges                = heapToDoubles(heapToDoubles(pCrackedEdges_mesh_out), bamgmeshout.CrackedEdgesSize);
+    bamgmeshout.VerticesOnGeomEdgeSize      = heapToint(pVerticesOnGeomEdgeSize_mesh_out, 2);
+    bamgmeshout.VerticesOnGeomEdge          = heapToDoubles(heapToDoubles(pVerticesOnGeomEdge_mesh_out), bamgmeshout.VerticesOnGeomEdgeSize);
+    bamgmeshout.VerticesOnGeomVertexSize    = heapToint(pVerticesOnGeomVertexSize_mesh_out, 2);
+    bamgmeshout.VerticesOnGeomVertex        = heapToDoubles(heapToDoubles(pVerticesOnGeomVertex_mesh_out), bamgmeshout.VerticesOnGeomVertexSize);
+    bamgmeshout.EdgesOnGeomEdgeSize         = heapToint(pEdgesOnGeomEdgeSize_mesh_out, 2);
+    bamgmeshout.EdgesOnGeomEdge             = heapToDoubles(heapToDoubles(pEdgesOnGeomEdge_mesh_out), bamgmeshout.EdgesOnGeomEdgeSize);
+    bamgmeshout.IssmSegmentsSize            = heapToint(pIssmSegmentsSize_mesh_out, 2);
+    bamgmeshout.IssmSegments                = heapToDoubles(heapToDoubles(pIssmSegments_mesh_out), bamgmeshout.IssmSegmentsSize);
+
+    bamggeomout.VerticesSize                = heapToint(pVerticesSize_mesh_out, 2);
+    bamggeomout.Vertices                    = heapToDoubles(heapToDoubles(pVertices_mesh_out), bamggeomout.VerticesSize);
+    bamggeomout.EdgesSize                   = heapToint(pEdgesSize_mesh_out, 2);
+    bamggeomout.Edges                       = heapToDoubles(heapToDoubles(pEdges_mesh_out), bamggeomout.EdgesSize);
+    bamggeomout.CornersSize                 = heapToint(pCornersSize_mesh_out, 2);
+    bamggeomout.Corners                     = heapToDoubles(heapToDoubles(pCorners_mesh_out), bamggeomout.CornersSize);
+    bamggeomout.RequiredVerticesSize        = heapToint(pRequiredVerticesSize_mesh_out, 2);
+    bamggeomout.RequiredVertices            = heapToDoubles(heapToDoubles(pRequiredVertices_mesh_out), bamggeomout.RequiredVerticesSize);
+    bamggeomout.RequiredCrackedEdgesSize    = heapToint(pRequiredCrackedEdgesSize_mesh_out, 2);
+    bamggeomout.RequiredCrackedEdges        = heapToDoubles(heapToDoubles(pRequiredCrackedEdges_mesh_out), bamggeomout.RequiredCrackedEdgesSize);
+    bamggeomout.CrackedEdgesSize            = heapToint(pCrackedEdgesSize_mesh_out, 2);
+    bamggeomout.CrackedEdges                = heapToDoubles(heapToDoubles(pCrackedEdges_mesh_out), bamggeomout.CrackedEdgesSize);
+    bamggeomout.SubDomainSize               = heapToint(pSubDomainsSize_mesh_out, 2);
+    bamggeomout.SubDomains                  = heapToDoubles(heapToDoubles(pSubDomains_mesh_out), bamggeomout.SubDomainsSize);
+	/*}}}*/
+
+	var return_array=[bamgmeshout, bamggeomout];
+
+	/*Free ressources: */
+    Module._free(pVertices_mesh_in);
+    Module._free(pVerticesSize_mesh_in);
+    Module._free(pEdges_mesh_in);
+    Module._free(pEdgesSize_mesh_in);
+    Module._free(pTriangles_mesh_in);
+    Module._free(pTrianglesSize_mesh_in);
+    Module._free(pCrackedEdges_mesh_in);
+    Module._free(pCrackedEdgesSize_mesh_in);
+    Module._free(pVerticesOnGeomEdge_mesh_in);
+    Module._free(pVerticesOnGeomEdgeSize_mesh_in);
+    Module._free(pVerticesOnGeomVertex_mesh_in);
+    Module._free(pVerticesOnGeomVertexSize_mesh_in);
+    Module._free(pEdgesOnGeomEdge_mesh_in);
+    Module._free(pEdgesOnGeomEdgeSize_mesh_in);
+    Module._free(pIssmSegments_mesh_in);
+    Module._free(pIssmSegmentsSize_mesh_in);
+
+    Module._free(pVertices_geom_in);
+    Module._free(pVerticesSize_geom_in);
+    Module._free(pEdges_geom_in);
+    Module._free(pEdgesSize_geom_in);
+    Module._free(pCorners_geom_in);
+    Module._free(pCornersSize_geom_in);
+    Module._free(pRequiredVertices_geom_in);
+    Module._free(pRequiredVerticesSize_geom_in);
+    Module._free(pRequiredCrackedEdges_geom_in);
+    Module._free(pRequiredCrackedEdgesSize_geom_in);
+    Module._free(pCrackedEdges_geom_in);
+    Module._free(pCrackedEdgesSize_geom_in);
+    Module._free(pSubDomains_geom_in);
+    Module._free(pSubDomainsSize_geom_in);
+
+    Module._free(pVertices_geom_out);
+    Module._free(pVerticesSize_geom_out);
+    Module._free(pEdges_geom_out);
+    Module._free(pEdgesSize_geom_out);
+    Module._free(pCorners_geom_out);
+    Module._free(pCornersSize_geom_out);
+    Module._free(pRequiredVertices_geom_out);
+    Module._free(pRequiredVerticesSize_geom_out);
+    Module._free(pRequiredCrackedEdges_geom_out);
+    Module._free(pRequiredCrackedEdgesSize_geom_out);
+    Module._free(pCrackedEdges_geom_out);
+    Module._free(pCrackedEdgesSize_geom_out);
+    Module._free(pSubDomains_geom_out);
+    Module._free(pSubDomainsSize_geom_out);
+
+    Module._free(pVertices_mesh_out);
+    Module._free(pVerticesSize_mesh_out);
+    Module._free(pEdges_mesh_out);
+    Module._free(pEdgesSize_mesh_out);
+    Module._free(pTriangles_mesh_out);
+    Module._free(pTrianglesSize_mesh_out);
+    Module._free(pCrackedEdges_mesh_out);
+    Module._free(pCrackedEdgesSize_mesh_out);
+    Module._free(pVerticesOnGeomEdge_mesh_out);
+    Module._free(pVerticesOnGeomEdgeSize_mesh_out);
+    Module._free(pVerticesOnGeomVertex_mesh_out);
+    Module._free(pVerticesOnGeomVertexSize_mesh_out);
+    Module._free(pEdgesOnGeomEdge_mesh_out);
+    Module._free(pEdgesOnGeomEdgeSize_mesh_out);
+    Module._free(pIssmSegments_mesh_out);
+    Module._free(pIssmSegmentsSize_mesh_out);
+
+    Module._free(phminVertices);
+    Module._free(phminVerticesSize);
+    Module._free(phmaxVertices);
+    Module._free(phmaxVerticesSize);
+    Module._free(phVertices);
+    Module._free(phVerticesSize);
+    Module._free(pmetric);
+    Module._free(pmetricSize);
+    Module._free(pfield);
+    Module._free(pfieldSize);
+
+	return return_array;
+}
Index: /issm/trunk-jpl/src/wrappers/javascript/io/FetchJavascriptData.cpp
===================================================================
--- /issm/trunk-jpl/src/wrappers/javascript/io/FetchJavascriptData.cpp	(revision 22876)
+++ /issm/trunk-jpl/src/wrappers/javascript/io/FetchJavascriptData.cpp	(revision 22877)
@@ -122,4 +122,82 @@
 	
 	*pcontours=contours;
+}
+/*}}}*/
+/*FUNCTION FetchData(BamgGeom** pbamggeom, double* vertices, int nods){{{*/
+void FetchData(BamgGeom** pbamggeom, double* Vertices, int* VerticesSize, double* Edges, int* EdgesSize, double* Corners, int* CornersSize, double* RequiredVertices, int* RequiredVerticesSize, double* RequiredCrackedEdges, int* RequiredCrackedEdgesSize, double* CrackedEdges, int* CrackedEdgesSize, double* SubDomains, int* SubDomainsSize){
+
+	/*Initialize output*/
+	BamgGeom* bamggeom = new BamgGeom();
+
+	/*Fetch all fields*/
+	FetchData(&bamggeom->Vertices, &bamggeom->VerticesSize[0], &bamggeom->VerticesSize[1], Vertices, VerticesSize[0], VerticesSize[1]);
+	FetchData(&bamggeom->Edges, &bamggeom->EdgesSize[0], &bamggeom->EdgesSize[1], Edges, EdgesSize[0], EdgesSize[1]);
+	FetchData(&bamggeom->Corners, &bamggeom->CornersSize[0], &bamggeom->CornersSize[1], Corners, CornersSize[0], CornersSize[1]);
+	FetchData(&bamggeom->RequiredVertices, &bamggeom->RequiredVerticesSize[0], &bamggeom->RequiredVerticesSize[1], RequiredVertices, RequiredVerticesSize[0], RequiredVerticesSize[1]);
+	FetchData(&bamggeom->RequiredEdges, &bamggeom->RequiredEdgesSize[0], &bamggeom->RequiredEdgesSize[1], RequiredEdges, RequiredEdgesSize[0], RequiredEdgesSize[1]);
+	FetchData(&bamggeom->CrackedEdges, &bamggeom->CrackedEdgesSize[0], &bamggeom->CrackedEdgesSize[1], CrackedEdges, CrackedEdgesSize[0], CrackedEdgesSize[1]);
+	FetchData(&bamggeom->SubDomains, &bamggeom->SubDomainsSize[0], &bamggeom->SubDomainsSize[1], SubDomains, SubDomainsSize[0], SubDomainsSize[1]);
+
+	/*Assign output pointers:*/
+	*pbamggeom = bamggeom;
+}
+/*}}}*/
+/*FUNCTION FetchData(BamgMesh** pbamgmesh, double* vertices, int nods){{{*/
+void FetchData(BamgMesh** pbamgmesh, double* Vertices, int* VerticesSize, double* Edges, int* EdgesSize, double* Triangles, int* TrianglesSize, double* CrackedEdges, int* CrackedEdgesSize, double* VerticesOnGeomEdge, int* VerticesOnGeomEdgeSize, double* VerticesOnGeomVertex, int* VerticesOnGeomVertexSize, double* EdgesOnGeomEdge, int* EdgesOnGeomEdgeSize, double* IssmSegments, int* IssmSegmentsSize){
+
+	/*Initialize output*/
+	BamgMesh* bamgmesh = new BamgMesh();
+
+	/*Fetch all fields*/
+	FetchData(&bamggeom->Vertices, &bamggeom->VerticesSize[0], &bamggeom->VerticesSize[1], Vertices, VerticesSize[0], VerticesSize[1]);
+	FetchData(&bamggeom->Edges, &bamggeom->EdgesSize[0], &bamggeom->EdgesSize[1], Edges, EdgesSize[0], EdgesSize[1]);
+	FetchData(&bamggeom->Triangles, &bamggeom->TrianglesSize[0], &bamggeom->TrianglesSize[1], Triangles, TrianglesSize[0], TrianglesSize[1]);
+	FetchData(&bamggeom->CrackedEdges, &bamggeom->CrackedEdgesSize[0], &bamggeom->CrackedEdgesSize[1], CrackedEdges, CrackedEdgesSize[0], CrackedEdgesSize[1]);
+	FetchData(&bamggeom->VerticesOnGeomEdge, &bamggeom->VerticesOnGeomEdgeSize[0], &bamggeom->VerticesOnGeomEdgeSize[1], VerticesOnGeomEdge, VerticesOnGeomEdgeSize[0], VerticesOnGeomEdgeSize[1]);
+	FetchData(&bamggeom->VerticesOnGeomVertex, &bamggeom->VerticesOnGeomVertexSize[0], &bamggeom->VerticesOnGeomVertexSize[1], VerticesOnGeomVertex, VerticesOnGeomVertexSize[0], VerticesOnGeomVertexSize[1]);
+	FetchData(&bamggeom->EdgesOnGeomEdge, &bamggeom->EdgesOnGeomEdgeSize[0], &bamggeom->EdgesOnGeomEdgeSize[1], EdgesOnGeomEdge, EdgesOnGeomEdgeSize[0], EdgesOnGeomEdgeSize[1]);
+	FetchData(&bamggeom->IssmSegments, &bamggeom->IssmSegmentsSize[0], &bamggeom->IssmSegmentsSize[1], IssmSegments, IssmSegmentsSize[0], IssmSegmentsSize[1]);
+
+	/*Assign output pointers:*/
+	*pbamgmesh = bamgmesh;
+}
+/*}}}*/
+/*FUNCTION FetchData(BamgMesh** pbamgopts, double* vertices, int nods){{{*/
+void FetchData(BamgOpts** pbamgopts, double anisomax, double cutoff, double coeff, double errg, double gradation, int Hessiantype, int maxnbv, double maxsubdiv, int Metrictype, int nbjacobi, int nbsmooth, double omega, double power, int verbose, int Crack, int KeepVertices, int splitcorners, double hmin, double hmax, double* hminVertices, double* hmaxVertices, double* hVertices, double* metric, double* field, double err){
+
+	BamgOpts *bamgopts      = new BamgOpts();
+
+    /*Parameters*/
+    bamgopts->anisomax	    = anisomax;
+    bamgopts->coeff	        = coeff;
+    bamgopts->cutoff    	= cutoff;
+    bamgopts->errg	        = errg;
+    bamgopts->gradation	    = gradation;
+    bamgopts->Hessiantype	= Hessiantype;
+    bamgopts->maxnbv	    = maxnbv;
+    bamgopts->maxsubdiv	    = maxsubdiv;
+    bamgopts->Metrictype	= Metrictype;
+    bamgopts->nbjacobi	    = nbjacobi;
+    bamgopts->nbsmooth	    = nbsmooth;
+    bamgopts->omega	        = omega;
+    bamgopts->power	        = power;
+    bamgopts->verbose	    = verbose;
+
+    /*Flags*/
+    bamgopts->Crack	        = Crack;
+    bamgopts->KeepVertices	= KeepVertices;
+    bamgopts->splitcorners	= splitcorners;
+        
+    /*Metric related*/
+    bamgopts->hmin	        = hmin;
+    bamgopts->hmax       	= hmax;
+    bamgopts->hminVertices	= xNew<IssmPDouble>(hminVerticesNods);  xMemCpy<IssmPDouble>(bamgopts->hminVertices, hminVertices, hminVerticesNods);
+    bamgopts->hmaxVertices  = xNew<IssmPDouble>(hmaxVerticesNods);  xMemCpy<IssmPDouble>(bamgopts->hmaxVertices, hmaxVertices, hmaxVerticesNods);
+    bamgopts->hVertices     = xNew<IssmPDouble>(hVerticesNods);     xMemCpy<IssmPDouble>(bamgopts->hVertices, hVertices, hVerticesNods);
+    bamgopts->field	        = xNew<IssmPDouble>(fieldNods);         xMemCpy<IssmPDouble>(bamgopts->field, field, fieldNods);
+    bamgopts->metric	    = xNew<IssmPDouble>(metricNods);        xMemCpy<IssmPDouble>(bamgopts->metric, metric, metricNods);
+    bamgopts->err	        = err;
+
+	*pbamgopts              = bamgopts;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/wrappers/javascript/io/javascriptio.h
===================================================================
--- /issm/trunk-jpl/src/wrappers/javascript/io/javascriptio.h	(revision 22876)
+++ /issm/trunk-jpl/src/wrappers/javascript/io/javascriptio.h	(revision 22877)
@@ -37,4 +37,7 @@
 void FetchData(int **pmatrix, int* pM, int* pN, int* matrixin, int M, int N);
 void FetchData(Contours** pcontours,double* x, double* y, int nods);
+void FetchData(BamgGeom** pbamggeom, double* Vertices, int* VerticesSize, double* Edges, int* EdgesSize, double* Corners, int* CornersSize, double* RequiredVertices, int* RequiredVerticesSize, double* RequiredCrackedEdges, int* RequiredCrackedEdgesSize, double* CrackedEdges, int* CrackedEdgesSize, double* SubDomains, int* SubDomainsSize);
+void FetchData(BamgMesh** pbamgmesh, double* Vertices, int* VerticesSize, double* Edges, int* EdgesSize, double* Triangles, int* TrianglesSize, double* CrackedEdges, int* CrackedEdgesSize, double* VerticesOnGeomEdge, int* VerticesOnGeomEdgeSize, double* VerticesOnGeomVertex, int* VerticesOnGeomVertexSize, double* EdgesOnGeomEdge, int* EdgesOnGeomEdgeSize, double* IssmSegments, int* IssmSegmentsSize);
+void FetchData(BamgOpts** pbamgopts, double anisomax, double cutoff, double coeff, double errg, double gradation, int Hessiantype, int maxnbv, double maxsubdiv, int Metrictype, int nbjacobi, int nbsmooth, double omega, double power, int verbose, int Crack, int KeepVertices, int splitcorners, double hmin, double hmax, double* hminVertices, double* hmaxVertices, double* hVertices, double* metric, double* field, double err);
 void FetchData(Options** poptions,int NRHS, int nrhs, const char* optionname, double optionvalue);
 void FetchData(int* pinteger,int integer);
