Index: /issm/trunk-jpl/src/m/Makefile.am =================================================================== --- /issm/trunk-jpl/src/m/Makefile.am (revision 22097) +++ /issm/trunk-jpl/src/m/Makefile.am (revision 22098) @@ -97,4 +97,5 @@ ${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 \ @@ -107,4 +108,5 @@ ${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 \ @@ -131,11 +133,8 @@ ${ISSM_DIR}/src/m/consistency/checkfield.js \ ${ISSM_DIR}/src/m/consistency/ismodelselfconsistent.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/io/fileptr.js \ - ${ISSM_DIR}/src/m/io/savemodel.js \ - ${ISSM_DIR}/src/m/io/loadmodel.js \ - ${ISSM_DIR}/src/m/io/writetofile.js \ - ${ISSM_DIR}/src/m/io/writetofile.js \ + ${ISSM_DIR}/src/m/io/*.js \ ${ISSM_DIR}/src/m/materials/paterson.js \ ${ISSM_DIR}/src/m/mesh/triangle.js \ @@ -145,18 +144,5 @@ ${ISSM_DIR}/src/m/parameterization/setflowequation.js \ ${ISSM_DIR}/src/m/parameterization/setmask.js \ - ${ISSM_DIR}/src/m/plot/slider.js \ - ${ISSM_DIR}/src/m/plot/applyoptions.js \ - ${ISSM_DIR}/src/m/plot/checkplotoptions.js \ - ${ISSM_DIR}/src/m/plot/plot_manager.js \ - ${ISSM_DIR}/src/m/plot/plot_mesh.js \ - ${ISSM_DIR}/src/m/plot/plot_overlay.js \ - ${ISSM_DIR}/src/m/plot/plot_unit.js \ - ${ISSM_DIR}/src/m/plot/plot_quiver.js \ - ${ISSM_DIR}/src/m/plot/plotdoc.js \ - ${ISSM_DIR}/src/m/plot/plotmodel.js \ - ${ISSM_DIR}/src/m/plot/processdata.js \ - ${ISSM_DIR}/src/m/plot/processmesh.js \ - ${ISSM_DIR}/src/m/plot/gauge.js \ - ${ISSM_DIR}/src/m/plot/webgl.js \ + ${ISSM_DIR}/src/m/plot/*.js \ ${ISSM_DIR}/src/m/print/sprintf.js \ ${ISSM_DIR}/src/m/solve/loadresultsfrombuffer.js \ Index: /issm/trunk-jpl/src/m/array/arrayoperations.js =================================================================== --- /issm/trunk-jpl/src/m/array/arrayoperations.js (revision 22097) +++ /issm/trunk-jpl/src/m/array/arrayoperations.js (revision 22098) @@ -93,4 +93,27 @@ 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) { //{{{ + // Make sure that both typed arrays are of the same type + if(Object.prototype.toString.call(a) !== Object.prototype.toString.call(b)) { + throw 'The types of the two arguments passed for parameters a and b do not match.'; + } + var c; + if(a instanceof Float64Array) { + c = new a.constructor(a.length + b.length); + c.set(a); + c.set(b, a.length); + } + else { + c = clone(a).concat(clone(b)); + } + return c; } //}}} function ListToMatrix(list, elementsPerSubArray) { //{{{ @@ -181,5 +204,5 @@ for(var i=0;i',0); + md = checkfield(md,'fieldname','love.frequencies','NaN',1,'Inf',1,'numel',md.love.nfreq); + md = checkfield(md,'fieldname','love.sh_nmax','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.sh_nmin','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.g0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.r0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',1,'>',0); + md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0, 1]); + md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12); + if (md.love.sh_nmin<=1 && md.love.forcing_type==9) { + throw 'Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.'; + } + } // }}} + this.marshall = function(md,prefix,fid) { // {{{ + + WriteData(fid,prefix,'object',this,'fieldname','nfreq','format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','frequencies','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',this,'fieldname','sh_nmax','format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','sh_nmin','format','Integer'); + WriteData(fid,prefix,'object',this,'fieldname','g0','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','r0','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','mu0','format','Double'); + WriteData(fid,prefix,'object',this,'fieldname','allow_layer_deletion','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','forcing_type','format','Integer'); + + } // }}} + //properties + // {{{ + + this.nfreq = NaN; + this.frequencies = NaN; + this.sh_nmax = NaN; + this.sh_nmin = NaN; + this.g0 = NaN; + this.r0 = NaN; + this.mu0 = NaN; + this.allow_layer_deletion = NaN; + this.forcing_type = NaN; + + //set defaults + this.setdefaultparameters(); + // }}} +} Index: /issm/trunk-jpl/src/m/classes/inversion.js =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.js (revision 22097) +++ /issm/trunk-jpl/src/m/classes/inversion.js (revision 22098) @@ -81,11 +81,11 @@ this.thickness_obs=project3d(md,'vector',this.thickness_obs,'type','node'); - if (numel(this.cost_functions_coefficients)>1) { + if (this.cost_functions_coefficients.length>1) { this.cost_functions_coefficients=project3d(md,'vector',this.cost_functions_coefficients,'type','node'); } - if (numel(this.min_parameters)>1) { + if (this.min_parameters.length>1) { this.min_parameters=project3d(md,'vector',this.min_parameters,'type','node'); } - if (numel(this.max_parameters)>1) { + if (this.max_parameters.length>1) { this.max_parameters=project3d(md,'vector',this.max_parameters,'type','node'); } Index: /issm/trunk-jpl/src/m/classes/model.js =================================================================== --- /issm/trunk-jpl/src/m/classes/model.js (revision 22097) +++ /issm/trunk-jpl/src/m/classes/model.js (revision 22098) @@ -34,8 +34,9 @@ console.log(sprintf("//19s: //-22s -- //s","thermal" ,"[1x1 " + typeof(this.thermal) + "]","parameters for thermal solution")); console.log(sprintf("//19s: //-22s -- //s","steadystate" ,"[1x1 " + typeof(this.steadystate) + "]","parameters for steadystate solution")); - console.log(sprintf("//19s: //-22s -- //s","trans" ,"[1x1 " + typeof(this.trans) + "]","parameters for trans solution")); + console.log(sprintf("//19s: //-22s -- //s","trans" ,"[1x1 " + typeof(this.trans) + "]","parameters for trans solution")); console.log(sprintf("//19s: //-22s -- //s","levelset" ,"[1x1 " + typeof(this.levelset) + "]","parameters for moving boundaries (level-set method)")); console.log(sprintf("//19s: //-22s -- //s","calving" ,"[1x1 " + typeof(this.calving) + "]","parameters for calving")); console.log(sprintf("//19s: //-22s -- //s","gia" ,"[1x1 " + typeof(this.gia) + "]","parameters for gia solution")); + console.log(sprintf("//19s: //-22s -- //s","love" ,"[1x1 " + typeof(this.love) + "]","parameters for love solution")); console.log(sprintf("//19s: //-22s -- //s","esa" ,"[1x1 " + typeof(this.esa) + "]","parameters for elastic adjustment solution")); console.log(sprintf("//19s: //-22s -- //s","autodiff" ,"[1x1 " + typeof(this.autodiff) + "]","automatic differentiation parameters")); @@ -55,15 +56,15 @@ this.constants = new constants(); this.geometry = new geometry(); - this.initialization = new initialization(); this.smb = new SMBforcing(); this.basalforcings = new basalforcings(); + this.materials = new matice(); + this.damage = new damage(); this.friction = new friction(); + this.flowequation = new flowequation(); + this.timestepping = new timestepping(); + this.initialization = new initialization(); this.rifts = new rifts(); this.slr = new slr(); - this.timestepping = new timestepping(); - this.groundingline = new groundingline(); - this.materials = new matice(); - this.damage = new damage(); - this.flowequation = new flowequation(); + this.debug = new debug(); this.verbose = new verbose(); @@ -71,6 +72,8 @@ this.toolkits = new toolkits(); this.cluster = new local(); + this.balancethickness = new balancethickness(); this.stressbalance = new stressbalance(); + this.groundingline = new groundingline(); this.hydrology = new hydrologyshreve(); this.masstransport = new masstransport(); @@ -80,5 +83,7 @@ this.levelset = new levelset(); this.calving = new calving(); + this.love = new fourierlove(); this.gia = new giaivins(); + this.esa = new esa(); this.autodiff = new autodiff(); @@ -86,7 +91,8 @@ this.qmu = new qmu(); this.amr = new amr(); - this.radaroverlay = new radaroverlay(); + this.results = {}; this.outputdefinition = new outputdefinition(); + this.radaroverlay = new radaroverlay(); this.miscellaneous = new miscellaneous(); this.priv = new priv(); @@ -113,467 +119,435 @@ } //}}} - this.extrude = function(md) { //{{{ - //EXTRUDE - vertically extrude a 2d mesh - // - // vertically extrude a 2d mesh and create corresponding 3d mesh. - // The vertical distribution can: - // - follow a polynomial law - // - follow two polynomial laws, one for the lower part and one for the upper part of the mesh - // - be discribed by a list of coefficients (between 0 and 1) - // - // - // Usage: - // md=extrude(md,numlayers,extrusionexponent); - // md=extrude(md,numlayers,lowerexponent,upperexponent); - // md=extrude(md,listofcoefficients); - // - // Example: - // md=extrude(md,15,1.3); - // md=extrude(md,15,1.3,1.2); - // md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]); - // - // See also: MODELEXTRACT, COLLAPSE - - //some checks on list of arguments - var argc = arguments.length; - - if ((argc > 4) | (argc < 2)) { - throw "extrude error message"; - } - - //Extrude the mesh - if (argc==2) { //list of coefficients - clist=arguments[0]; - - if (clist.some(function(x) { return x < 0; }) | clist.some(function(x) { return x > 1; })) - throw 'extrusioncoefficients must be between 0 and 1'; - clistcopy.push.apply(clistcopy, [0, 1]); - - extrusionlist=clistcopy.reverse().filter(function(e, i, arr) { //only keep unique elements - return arr.indexOf(e, i+1) === -1; - }).sort(); - - numlayers=extrusionlist.length; - } else if (argc==3) { //one polynomial law - if (arguments[1]<=0) { - //help extrude; - throw 'extrusionexponent must be >=0'; - } - numlayers=arguments[0]; - - extrusionlist = []; - - for (var i = 0; i < numlayers; i++) { - extrusionlist.push(Math.pow(i/(numlayers-1), arguments[1])); - }; - - } else if (argc==4) { //two polynomial laws - numlayers=arguments[0]; - lowerexp=arguments[1]; - upperexp=arguments[2]; - - if (arguments[1]<=0 | args[2]<=0) { - throw 'lower and upper extrusionexponents must be >=0'; - } - - lowerextrusionlist = []; - for (var i = 0; i <= 1; i += (2/(numlayers-1))) { - lowerextrusionlist.push(Math.pow(i, lowerexp)/2) - }; - - upperextrusionlist = []; - for (var i = 0; i <= 1; i += (2/(numlayers-1))) { - upperextrusionlist.push(Math.pow(i, upperexp)/2) - }; - - //lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2; - //upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2; - - var temp = lowerextrusionlist; - temp.push.apply(temp, [1]); - - var temp2 = upperextrusionlist.map(function(x) { return -x; }); - temp.push.apply(temp, temp2); - - clistcopy.push.apply(clistcopy, [0, 1]); - //extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist])); - extrusionlist=temp.sort(); - } - - if (numlayers<2) { - console.error('number of layers should be at least 2'); - } - - if (md.mesh.domaintype() === '3D') { - console.error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); - } - - //Initialize with the 2d mesh - var mesh2d = md.mesh; - md.mesh=new mesh3dprisms(); - - md.mesh.x = mesh2d.x; - md.mesh.y = mesh2d.y; - md.mesh.elements = mesh2d.elements; - md.mesh.numberofelements = mesh2d.numberofelements; - md.mesh.numberofvertices = mesh2d.numberofvertices; - - md.mesh.lat = mesh2d.lat; - md.mesh.long = mesh2d.long; - md.mesh.epsg = mesh2d.epsg; - - md.mesh.vertexonboundary = mesh2d.vertexonboundary; - md.mesh.vertexconnectivity = mesh2d.vertexconnectivity; - md.mesh.elementconnectivity = mesh2d.elementconnectivity; - md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity; - - md.mesh.extractedvertices = mesh2d.extractedvertices; - md.mesh.extractedelements = mesh2d.extractedelements; - - x3d=[]; - y3d=[]; - z3d=[]; //the lower node is on the bed - thickness3d=md.geometry.thickness; //thickness and bed for these nodes - bed3d=md.geometry.base; - - //Create the new layers - for (var i = 1; i <= numlayers; i++) { - x3d=[x3d, md.mesh.x]; - y3d=[y3d, md.mesh.y]; - //nodes are distributed between bed and surface accordingly to the given exponent - z3d=[z3d, bed3d+thickness3d*extrusionlist(i)]; - }; - - number_nodes3d=x3d.length; //number of 3d nodes for the non extruded part of the mesh - - //Extrude elements - elements3d=[]; - for (var i = 0; i < numlayers; i++) { - elements3d.push(md.mesh.elements+(i-1)*md.mesh.numberofvertices, md.mesh.elements+i*md.mesh.numberofvertices); //Create the elements of the 3d mesh for the non extruded part - }; - - number_el3d=elements3d.length; //number of 3d nodes for the non extruded part of the mesh - - function fillArray(value, len) { - if (len == 0) return []; - var a = [value]; - while (a.length * 2 <= len) a = a.concat(a); - if (a.length < len) a = a.concat(a.slice(0, len - a.length)); - return a; - } - - //Keep a trace of lower and upper nodes - //lowervertex=NaN*ones(number_nodes3d,1); - //uppervertex=NaN*ones(number_nodes3d,1); - - lowervertex = fillArray(NaN, number_nodes3d); - uppervertex = fillArray(NaN, number_nodes3d); - - //lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices; - //uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d; - - for (var i = md.mesh.numberofvertices+1, k = 1; i < lowervertex.length && k <= (numlayers-1)*md.mesh.numberofvertices; i++, k++) { - lowervertex[i] = k; - }; - - for (var i = 1, k = md.mesh.numberofvertices+1; i <= (numlayers-1)*md.mesh.numberofvertices && k <= number_nodes3d; i++, k++) { - uppervertex[i] = k; - }; - - md.mesh.lowervertex=lowervertex; - md.mesh.uppervertex=uppervertex; - - //same for lower and upper elements - lowerelements = fillArray(NaN, number_el3d); - upperelements = fillArray(NaN, number_el3d); - //lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements; - //upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements; - - for (var i = md.mesh.numberofvertices+1, k = 1; i < lowerelements.length && k <= (numlayers-2)*md.mesh.numberofelements; i++, k++) { - lowerelements[i] = k; - }; - - for (var i = 1, k = md.mesh.numberofelements + 1; i <= (numlayers-2)*md.mesh.numberofelements && k <= (numlayers-1)*md.mesh.numberofelements; i++, k++) { - upperelements[i] = k; - }; - - md.mesh.lowerelements=lowerelements; - md.mesh.upperelements=upperelements; - - //Save old mesh - md.mesh.x2d=md.mesh.x; - md.mesh.y2d=md.mesh.y; - md.mesh.elements2d=md.mesh.elements; - md.mesh.numberofelements2d=md.mesh.numberofelements; - md.mesh.numberofvertices2d=md.mesh.numberofvertices; - - //Build global 3d mesh - md.mesh.elements=elements3d; - md.mesh.x=x3d; - md.mesh.y=y3d; - md.mesh.z=z3d; - md.mesh.numberofelements=number_el3d; - md.mesh.numberofvertices=number_nodes3d; - md.mesh.numberoflayers=numlayers; - - //Ok, now deal with the other fields from the 2d mesh: - - //bedinfo and surface info - md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1); - md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); - md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); - - //lat long - md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node'); - md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node'); - - //md.geometry=md.geometry.extrude(md); - //md.friction=md.friction.extrude(md); - //md.inversion=md.inversion.extrude(md); - //md.smb=md.smb.extrude(md); - //md.initialization=md.initialization.extrude(md); - - //md.flowequation=md.flowequation.extrude(md); - //md.stressbalance=md.stressbalance.extrude(md); - //md.thermal=md.thermal.extrude(md); - //md.masstransport=md.masstransport.extrude(md); - //md.levelset=md.levelset.extrude(md); - //md.calving=md.calving.extrude(md); - //md.hydrology = md.hydrology.extrude(md); - - ////connectivity - ////if ~isnan(md.mesh.elementconnectivity) - - //if (md.mesh.elementconnectivity.every(function(e, i, arr) { return e !== null && typeof e !== 'undefined'; })) { - //var temparr = []; - //for (var i = 0; i < numlayers; i++) { - //temparr.push(md.mesh.elementconnectivity); - //}; - - ////md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); - - //md.mesh.elementconnectivity = temparr; - - ////md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; - - //for (var i = 0; i < md.mesh.elementconnectivity.length; i++) { - //if (md.mesh.elementconnectivity[i] == 0) { - //md.mesh.elementconnectivity[i] = NaN; - //} - //}; - - //for (var i = 2; i < numlayers; i++) { - //for (var j = (i-1)*md.mesh.nuberofelements2d+1; i <= (i)*md.mesh.numberofelements2d; i++) { - //md.mesh.elementconnectivity[i] += md.mesh.numberofelements2d; - //}; - ////md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... - ////=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; - //}; - - //md.mesh.elementconnectivity.map(function(x) { return (isNaN(x)) ? 0 : x; }); - //} - - //md.materials=md.materials.extrude(md); - //md.damage=md.damage.extrude(md); - //md.mask=md.mask.extrude(md); - //md.qmu=md.qmu.extrude(md); - //md.basalforcings=md.basalforcings.extrude(md); - - ////increase connectivity if less than 25: - //if (md.mesh.average_vertex_connectivity<=25) - //md.mesh.average_vertex_connectivity=100; - } //}}} - this.collapse = function(md) { //{{{ - /* - *COLLAPSE - collapses a 3d mesh into a 2d mesh - * - * This routine collapses a 3d model into a 2d model - * and collapses all the fileds of the 3d model by - * taking their depth-averaged values - * - * Usage: - * md=collapse(md) - * - * See also: EXTRUDE, MODELEXTRACT - */ - - // Check that the model is really a 3d model - if (md.mesh.elementtype() !== 'Penta') { - console.error('collapse error message: only 3d mesh can be collapsed') - } - - // Start with changing all the fields from the 3d mesh - - // dealing with the friction law - // drag is limited to nodes that are on the bedrock. - if (md.friction.classname() === 'friction') { - md.friction.coefficient=project2d(md,md.friction.coefficient,1); - md.friction.p=project2d(md,md.friction.p,1); - md.friction.q=project2d(md,md.friction.q,1); - } else if (md.friction.classname() === 'frictionhydro') { - md.friction.q=project2d(md,md.friction.q,1); - md.friction.C=project2d(md,md.friction.C,1); - md.friction.As=project2d(md,md.friction.As,1); - md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1); - } else if (md.friction.classname() === 'frictionwaterlayer') { - md.friction.coefficient=project2d(md,md.friction.coefficient,1); - md.friction.p=project2d(md,md.friction.p,1); - md.friction.q=project2d(md,md.friction.q,1); - md.friction.water_layer=project2d(md,md.friction.water_layer,1); - } else if (md.friction.classname() === 'frictionweertman') { - md.friction.C=project2d(md,md.friction.C,1); - md.friction.m=project2d(md,md.friction.m,1); - } else if (md.friction.classname() === 'frictionweertmantemp') { - md.friction.C=project2d(md,md.friction.C,1); - md.friction.m=project2d(md,md.friction.m,1); - } else { - disp('friction type not supported'); - } - - // observations - if (!isNaN(md.inversion.vx_obs)) - md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); - - if (!isNaN(md.inversion.vy_obs)) - md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); - - if (!isNaN(md.inversion.vel_obs)) - md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); - - if (!isNaN(md.inversion.cost_functions_coefficients)) - md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); - - if (numel(md.inversion.min_parameters)>1) - md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); - - if (numel(md.inversion.max_parameters)>1) - md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); - - if (md.smb.classname() === 'SMBforcing' && !isNaN(md.smb.mass_balance)) { - md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers); - } else if (md.smb.classname() === 'SMBhenning' && !isNaN(md.smb.smbref)) { - md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers); - } - - // results - if (!isNaN(md.initialization.vx)) - md.initialization.vx=DepthAverage(md,md.initialization.vx); - if (!isNaN(md.initialization.vy)) - md.initialization.vy=DepthAverage(md,md.initialization.vy); - if (!isNaN(md.initialization.vz)) - md.initialization.vz=DepthAverage(md,md.initialization.vz); - if (!isNaN(md.initialization.vel)) - md.initialization.vel=DepthAverage(md,md.initialization.vel); - if (!isNaN(md.initialization.temperature)) - md.initialization.temperature=DepthAverage(md,md.initialization.temperature); - if (!isNaN(md.initialization.pressure)) - md.initialization.pressure=project2d(md,md.initialization.pressure,1); - if (!isNaN(md.initialization.sediment_head)) - md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1); - if (!isNaN(md.initialization.epl_head)) - md.initialization.epl_head=project2d(md,md.initialization.epl_head,1); - if (!isNaN(md.initialization.epl_thickness)) - md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1); - - // giaivins - if (!isNaN(md.gia.mantle_viscosity)) - md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); - if (!isNaN(md.gia.lithosphere_thickness)) - md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1); - - // elementstype - if (!isNaN(md.flowequation.element_equation)) { - md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); - md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); - md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1); - md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1); - md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1); - } - - // boundary conditions - md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers); - md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers); - md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers); - md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers); - md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers); - md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers); - if (!isNaN(md.damage.spcdamage)) - md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); - md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); - - // Hydrologydc variables - if (md.hydrology.classname() === 'hydrologydc') { - md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1); - md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1); - md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1); - md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1); - - if(md.hydrology.isefficientlayer==1) - md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1); - } - - // materials - md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); - md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); - - // damage: - if (md.damage.isdamage) - md.damage.D=DepthAverage(md,md.damage.D); - - // special for thermal modeling: - if (!isNaN(md.basalforcings.groundedice_melting_rate)) - md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1); - - if (!isNaN(md.basalforcings.floatingice_melting_rate)) - md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1); - - md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); // bedrock only gets geothermal flux - - // update of connectivity matrix - md.mesh.average_vertex_connectivity=25; - - // Collapse the mesh - var nodes2d=md.mesh.numberofvertices2d; - var elements2d=md.mesh.numberofelements2d; - - // parameters - md.geometry.surface=project2d(md,md.geometry.surface,1); - md.geometry.thickness=project2d(md,md.geometry.thickness,1); - md.geometry.base=project2d(md,md.geometry.base,1); - if (!isNaN(md.geometry.bed)) - md.geometry.bed=project2d(md,md.geometry.bed,1); - - if (!isNaN(md.mask.groundedice_levelset)) - md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1); - - if (!isNaN(md.mask.ice_levelset)) - md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1); - - // lat long - if (numel(md.mesh.lat) === md.mesh.numberofvertices) - md.mesh.lat=project2d(md,md.mesh.lat,1); - if (numel(md.mesh.long) === md.mesh.numberofvertices) - md.mesh.long=project2d(md,md.mesh.long,1); - - // Initialize with the 2d mesh - var mesh = new mesh2d(); - mesh.x=md.mesh.x2d; - mesh.y=md.mesh.y2d; - mesh.numberofvertices=md.mesh.numberofvertices2d; - mesh.numberofelements=md.mesh.numberofelements2d; - mesh.elements=md.mesh.elements2d; - - if (!isNaN(md.mesh.vertexonboundary)) - mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); - if (!isNaN(md.mesh.elementconnectivity)) - mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); - - md.mesh=mesh; - md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); - md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); - md.mesh.segments=contourenvelope(md.mesh); - - return md; - } /*}}}*/ + this.extrude = function(md) { //{{{ + //EXTRUDE - vertically extrude a 2d mesh + // + // vertically extrude a 2d mesh and create corresponding 3d mesh. + // The vertical distribution can: + // - follow a polynomial law + // - follow two polynomial laws, one for the lower part and one for the upper part of the mesh + // - be described by a list of coefficients (between 0 and 1) + // + // + // Usage: + // md=extrude(md,numlayers,extrusionexponent); + // md=extrude(md,numlayers,lowerexponent,upperexponent); + // md=extrude(md,listofcoefficients); + // + // Example: + // md=extrude(md,15,1.3); + // md=extrude(md,15,1.3,1.2); + // md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]); + // + // See also: MODELEXTRACT, COLLAPSE + + //some checks on list of arguments + var argc = arguments.length; + var extrusionlist; + var numlayers; + + if ((argc > 4) | (argc < 2)) { + throw "extrude error message"; + } + + //Extrude the mesh + if (argc==2) { //list of coefficients + clist=arguments[0]; + if (ArrayAnyBelowStrict(clist,0) || ArrayAnyAboveStrict(clist,1)) { + throw 'extrusioncoefficients must be between 0 and 1'; + } + extrusionlist=ArraySort(ArrayUnique(clist.push(0,1))); + numlayers=extrusionlist.length; + } else if (argc==3) { //one polynomial law + if (arguments[1]<=0) { + throw 'extrusionexponent must be >=0'; + } + numlayers=arguments[1]; + extrusionlist = []; + for (var i = 0; i < numlayers; i++) { + extrusionlist.push(Math.pow(i/(numlayers-1), arguments[2])); + } + } else if (argc==4) { //two polynomial laws + numlayers=arguments[1]; + var lowerexp=arguments[2]; + var upperexp=arguments[3]; + + if (arguments[2]<=0 || arguments[3]<=0) { + throw 'lower and upper extrusionexponents must be >=0'; + } + + var lowerextrusionlist = []; + for (var i = 0; i <= 1; i += (2/(numlayers-1))) { + lowerextrusionlist.push(Math.pow(i, lowerexp)/2) + } + var upperextrusionlist = []; + for (var i = 0; i <= 1; i += (2/(numlayers-1))) { + upperextrusionlist.push(1-Math.pow(i, upperexp)/2) + } + extrusionlist=ArrayConcat(lowerextrusionlist,upperextrusionlist); + extrusionlist=ArraySort(ArrayUnique(extrusionlist.push(1))); + } + + if (numlayers<2) { + console.error('number of layers should be at least 2'); + } + if (md.mesh.domaintype() === '3D') { + console.error('Cannot extrude a 3d mesh (extrude cannot be called more than once)'); + } + + //Initialize with the 2d mesh + var mesh2d = md.mesh; + md.mesh=new mesh3dprisms(); + md.mesh.x = mesh2d.x; + md.mesh.y = mesh2d.y; + md.mesh.elements = mesh2d.elements; + md.mesh.numberofelements = mesh2d.numberofelements; + md.mesh.numberofvertices = mesh2d.numberofvertices; + + md.mesh.lat = mesh2d.lat; + md.mesh.long = mesh2d.long; + md.mesh.epsg = mesh2d.epsg; + + md.mesh.vertexonboundary = mesh2d.vertexonboundary; + md.mesh.vertexconnectivity = mesh2d.vertexconnectivity; + md.mesh.elementconnectivity = mesh2d.elementconnectivity; + md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity; + + md.mesh.extractedvertices = mesh2d.extractedvertices; + md.mesh.extractedelements = mesh2d.extractedelements; + + var x3d=new Float64Array(); + var y3d=new Float64Array(); + var z3d=new Float64Array(); //the lower node is on the bed + var thickness3d=md.geometry.thickness; //thickness and bed for these nodes + var bed3d=md.geometry.base; + + //Create the new layers + //Dimensions of x/y/z3d: md.mesh.numberofvertices * numlayers, 1 + for (var i = 1; i <= numlayers; i++) { + x3d=ArrayConcat(x3d, md.mesh.x); + y3d=ArrayConcat(y3d, md.mesh.y); + //nodes are distributed between bed and surface accordingly to the given exponent + z3d=ArrayConcat(z3d, FloatFix(thickness3d.map(function(value,index) { return bed3d[index] + value * extrusionlist[i-1]; }), thickness3d.length)); + } + var number_nodes3d=x3d.length; //number of 3d nodes for the non extruded part of the mesh + + //Extrude elements + //Create the elements of the 3d mesh for the non extruded part + //Dimensions of elements3d: md.mesh.numberofelements * (numlayers - 1), 6 + var elements3d=[]; + var elements_prisms=[]; + for (var i = 1; i < numlayers; i++) { + for (var j = 0; j < md.mesh.numberofelements; j++) { + elements_prisms = []; + for (var k = 0, offset = (i - 1) * md.mesh.numberofvertices; k < 3; k++) { + elements_prisms.push(md.mesh.elements[j][k]+offset); + } + for (var k = 0, offset = i * md.mesh.numberofvertices; k < 3; k++) { + elements_prisms.push(md.mesh.elements[j][k]+offset); + } + elements3d.push(elements_prisms); + } + } + number_el3d=elements3d.length; //number of 3d nodes for the non extruded part of the mesh + + //Keep a trace of lower and upper nodes + var lowervertex = NewArrayFill(number_nodes3d, NaN); + var uppervertex = NewArrayFill(number_nodes3d, NaN); + //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofvertices+1) + for (var i = md.mesh.numberofvertices, k = 1; i < lowervertex.length; i++, k++) { + lowervertex[i] = k; + }; + //Set last layer to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1) + for (var i = 0, k = md.mesh.numberofvertices+1; i < (numlayers-1)*md.mesh.numberofvertices; i++, k++) { + uppervertex[i] = k; + }; + md.mesh.lowervertex=lowervertex; + md.mesh.uppervertex=uppervertex; + + //same for lower and upper elements + var lowerelements = NewArrayFill(number_el3d, NaN); + var upperelements = NewArrayFill(number_el3d, NaN); + //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofelements+1) + for (var i = md.mesh.numberofelements, k = 1; i < lowerelements.length; i++, k++) { + lowerelements[i] = k; + }; + //Set last 2 layers to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1) + for (var i = 0, k = md.mesh.numberofelements + 1; i < (numlayers-2)*md.mesh.numberofelements; i++, k++) { + upperelements[i] = k; + }; + md.mesh.lowerelements=lowerelements; + md.mesh.upperelements=upperelements; + + //Save old mesh + md.mesh.x2d=md.mesh.x; + md.mesh.y2d=md.mesh.y; + md.mesh.elements2d=md.mesh.elements; + md.mesh.numberofelements2d=md.mesh.numberofelements; + md.mesh.numberofvertices2d=md.mesh.numberofvertices; + + //Build global 3d mesh + md.mesh.elements=elements3d; + md.mesh.x=x3d; + md.mesh.y=y3d; + md.mesh.z=z3d; + md.mesh.numberofelements=number_el3d; + md.mesh.numberofvertices=number_nodes3d; + md.mesh.numberoflayers=numlayers; + + //Ok, now deal with the other fields from the 2d mesh: + + //bedinfo and surface info + md.mesh.vertexonbase=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',1); + md.mesh.vertexonsurface=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers); + md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node'); + + //lat long + md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node'); + md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node'); + + md.geometry=md.geometry.extrude(md); + md.friction=md.friction.extrude(md); + md.inversion=md.inversion.extrude(md); + md.smb=md.smb.extrude(md); + md.initialization=md.initialization.extrude(md); + + md.flowequation=md.flowequation.extrude(md); + md.stressbalance=md.stressbalance.extrude(md); + md.thermal=md.thermal.extrude(md); + md.masstransport=md.masstransport.extrude(md); + md.levelset=md.levelset.extrude(md); + md.calving=md.calving.extrude(md); + md.hydrology = md.hydrology.extrude(md); + + //connectivity + if (!ArrayAnyNaN(md.mesh.elementconnectivity)) { + //md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); //Replicate the matrix across numlayers-1 + var temparr = []; + for (var i = 0; i < numlayers - 1; i++) { + temparr = ArrayConcat(temparr, md.mesh.elementconnectivity); + } + md.mesh.elementconnectivity = temparr; + + //md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN; + var indices = ArrayFind(md.mesh.elementconnectivity, 0); + for (var i = 0; i < indices.length; i++) { + md.mesh.elementconnectivity[i] = NaN; + }; + + for (var i = 2; i < numlayers; i++) { + //md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)... + //=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d; + for (var j = (i-1)*md.mesh.numberofelements2d; j <= (i)*md.mesh.numberofelements2d-1; j++) { + for (var k = 0; k < 3; k++) { + md.mesh.elementconnectivity[j][k] += md.mesh.numberofelements2d; + } + } + } + md.mesh.elementconnectivity = md.mesh.elementconnectivity.map(function(value) { return (Number.isNaN(value)) ? 0 : value; }); + } + + md.materials=md.materials.extrude(md); + md.damage=md.damage.extrude(md); + md.mask=md.mask.extrude(md); + md.qmu=md.qmu.extrude(md); + md.basalforcings=md.basalforcings.extrude(md); + + //increase connectivity if less than 25: + if (md.mesh.average_vertex_connectivity<=25) { + md.mesh.average_vertex_connectivity=100; + } + } //}}} + this.collapse = function(md) { //{{{ + /* + *COLLAPSE - collapses a 3d mesh into a 2d mesh + * + * This routine collapses a 3d model into a 2d model + * and collapses all the fileds of the 3d model by + * taking their depth-averaged values + * + * Usage: + * md=collapse(md) + * + * See also: EXTRUDE, MODELEXTRACT + */ + + // Check that the model is really a 3d model + if (md.mesh.elementtype() !== 'Penta') { + console.error('collapse error message: only 3d mesh can be collapsed') + } + + // Start with changing all the fields from the 3d mesh + + // dealing with the friction law + // drag is limited to nodes that are on the bedrock. + if (md.friction.classname() === 'friction') { + md.friction.coefficient=project2d(md,md.friction.coefficient,1); + md.friction.p=project2d(md,md.friction.p,1); + md.friction.q=project2d(md,md.friction.q,1); + } else if (md.friction.classname() === 'frictionhydro') { + md.friction.q=project2d(md,md.friction.q,1); + md.friction.C=project2d(md,md.friction.C,1); + md.friction.As=project2d(md,md.friction.As,1); + md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1); + } else if (md.friction.classname() === 'frictionwaterlayer') { + md.friction.coefficient=project2d(md,md.friction.coefficient,1); + md.friction.p=project2d(md,md.friction.p,1); + md.friction.q=project2d(md,md.friction.q,1); + md.friction.water_layer=project2d(md,md.friction.water_layer,1); + } else if (md.friction.classname() === 'frictionweertman') { + md.friction.C=project2d(md,md.friction.C,1); + md.friction.m=project2d(md,md.friction.m,1); + } else if (md.friction.classname() === 'frictionweertmantemp') { + md.friction.C=project2d(md,md.friction.C,1); + md.friction.m=project2d(md,md.friction.m,1); + } else { + disp('friction type not supported'); + } + + // observations + if (!Number.isNaN(md.inversion.vx_obs)) + md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); + + if (!Number.isNaN(md.inversion.vy_obs)) + md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); + + if (!Number.isNaN(md.inversion.vel_obs)) + md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); + + if (!Number.isNaN(md.inversion.cost_functions_coefficients)) + md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); + + if (numel(md.inversion.min_parameters)>1) + md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); + + if (numel(md.inversion.max_parameters)>1) + md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); + + if (md.smb.classname() === 'SMBforcing' && !Number.isNaN(md.smb.mass_balance)) { + md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers); + } else if (md.smb.classname() === 'SMBhenning' && !Number.isNaN(md.smb.smbref)) { + md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers); + } + + // results + if (!Number.isNaN(md.initialization.vx)) + md.initialization.vx=DepthAverage(md,md.initialization.vx); + if (!Number.isNaN(md.initialization.vy)) + md.initialization.vy=DepthAverage(md,md.initialization.vy); + if (!Number.isNaN(md.initialization.vz)) + md.initialization.vz=DepthAverage(md,md.initialization.vz); + if (!Number.isNaN(md.initialization.vel)) + md.initialization.vel=DepthAverage(md,md.initialization.vel); + if (!Number.isNaN(md.initialization.temperature)) + md.initialization.temperature=DepthAverage(md,md.initialization.temperature); + if (!Number.isNaN(md.initialization.pressure)) + md.initialization.pressure=project2d(md,md.initialization.pressure,1); + if (!Number.isNaN(md.initialization.sediment_head)) + md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1); + if (!Number.isNaN(md.initialization.epl_head)) + md.initialization.epl_head=project2d(md,md.initialization.epl_head,1); + if (!Number.isNaN(md.initialization.epl_thickness)) + md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1); + + // giaivins + if (!Number.isNaN(md.gia.mantle_viscosity)) + md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1); + if (!Number.isNaN(md.gia.lithosphere_thickness)) + md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1); + + // elementstype + if (!Number.isNaN(md.flowequation.element_equation)) { + md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1); + md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1); + md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1); + md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1); + md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1); + } + + // boundary conditions + md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers); + md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers); + md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers); + md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers); + md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers); + md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers); + if (!Number.isNaN(md.damage.spcdamage)) + md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers); + md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); + + // Hydrologydc variables + if (md.hydrology.classname() === 'hydrologydc') { + md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1); + md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1); + md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1); + md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1); + + if(md.hydrology.isefficientlayer==1) + md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1); + } + + // materials + md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); + md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); + + // damage: + if (md.damage.isdamage) + md.damage.D=DepthAverage(md,md.damage.D); + + // special for thermal modeling: + if (!Number.isNaN(md.basalforcings.groundedice_melting_rate)) + md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1); + + if (!Number.isNaN(md.basalforcings.floatingice_melting_rate)) + md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1); + + md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); // bedrock only gets geothermal flux + + // update of connectivity matrix + md.mesh.average_vertex_connectivity=25; + + // Collapse the mesh + var nodes2d=md.mesh.numberofvertices2d; + var elements2d=md.mesh.numberofelements2d; + + // parameters + md.geometry.surface=project2d(md,md.geometry.surface,1); + md.geometry.thickness=project2d(md,md.geometry.thickness,1); + md.geometry.base=project2d(md,md.geometry.base,1); + if (!Number.isNaN(md.geometry.bed)) + md.geometry.bed=project2d(md,md.geometry.bed,1); + + if (!Number.isNaN(md.mask.groundedice_levelset)) + md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1); + + if (!Number.isNaN(md.mask.ice_levelset)) + md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1); + + // lat long + if (numel(md.mesh.lat) === md.mesh.numberofvertices) + md.mesh.lat=project2d(md,md.mesh.lat,1); + if (numel(md.mesh.long) === md.mesh.numberofvertices) + md.mesh.long=project2d(md,md.mesh.long,1); + + // Initialize with the 2d mesh + var mesh = new mesh2d(); + mesh.x=md.mesh.x2d; + mesh.y=md.mesh.y2d; + mesh.numberofvertices=md.mesh.numberofvertices2d; + mesh.numberofelements=md.mesh.numberofelements2d; + mesh.elements=md.mesh.elements2d; + + if (!Number.isNaN(md.mesh.vertexonboundary)) + mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1); + if (!Number.isNaN(md.mesh.elementconnectivity)) + mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1); + + md.mesh=mesh; + md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); + md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); + md.mesh.segments=contourenvelope(md.mesh); + + return md; + } /*}}}*/ //properties // {{{ @@ -611,4 +585,5 @@ this.levelset = 0; this.calving = 0; + this.love = 0; this.gia = 0; this.esa = 0; Index: /issm/trunk-jpl/src/m/classes/qmu.js =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu.js (revision 22097) +++ /issm/trunk-jpl/src/m/classes/qmu.js (revision 22098) @@ -91,5 +91,5 @@ }// }}} this.extrude = function(md) {//{{{ - this.partition=project3d(md,'vector',math.transpose(this.partition),'type','node'); + if (!isNaN(this.partition)) this.partition=project3d(md,'vector',ArrayTranspose(this.partition),'type','node'); return this; }//}}} Index: /issm/trunk-jpl/src/m/classes/thermal.js =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.js (revision 22097) +++ /issm/trunk-jpl/src/m/classes/thermal.js (revision 22098) @@ -55,9 +55,26 @@ this.spctemperature=project3d(md,'vector',this.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN); if (md.initialization.temperature.length===md.mesh.numberofvertices) { - this.spctemperature = ArrayFill2D(Create2DArray(md.mesh.numberofvertices, 1), NaN); - var pos=ArrayFindNot2D(md.mesh.vertexonsurface, 0); - pos.forEach(function(x) { // impose observed temperature on surface - this.spctemperature[x] = md.initialization.temperature[x]; - }); + this.spctemperature = NewArrayFill(md.mesh.numberofvertices, NaN); + var pos=ArrayFindNot(md.mesh.vertexonsurface, 0); + // impose observed temperature on surface + for (var i=0,posIndex=0,count=0;i