Changeset 22098


Ignore:
Timestamp:
09/17/17 23:19:31 (8 years ago)
Author:
dlcheng
Message:

BUG: Fixing issues with test102.js extrusion functions. NEW: Adding fourierlove.js. CHG: Updating m/Makefile.am.

Location:
issm/trunk-jpl/src/m
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/Makefile.am

    r21945 r22098  
    9797                                ${ISSM_DIR}/src/m/classes/geometry.js \
    9898                                ${ISSM_DIR}/src/m/classes/giaivins.js \
     99                                ${ISSM_DIR}/src/m/classes/fourierlove.js \
    99100                                ${ISSM_DIR}/src/m/classes/esa.js \
    100101                                ${ISSM_DIR}/src/m/classes/groundingline.js \
    … …  
    107108                                ${ISSM_DIR}/src/m/classes/matice.js \
    108109                                ${ISSM_DIR}/src/m/classes/mesh2d.js \
     110                                ${ISSM_DIR}/src/m/classes/mesh3dprisms.js \
    109111                                ${ISSM_DIR}/src/m/classes/mesh3dsurface.js \
    110112                                ${ISSM_DIR}/src/m/classes/miscellaneous.js \
    … …  
    131133                                ${ISSM_DIR}/src/m/consistency/checkfield.js \
    132134                                ${ISSM_DIR}/src/m/consistency/ismodelselfconsistent.js \
     135                                ${ISSM_DIR}/src/m/extrusion/*.js \
    133136                                ${ISSM_DIR}/src/m/geometry/FlagElements.js \
    134137                                ${ISSM_DIR}/src/m/inversions/marshallcostfunctions.js \
    135                                 ${ISSM_DIR}/src/m/io/fileptr.js \
    136                                 ${ISSM_DIR}/src/m/io/savemodel.js \
    137                                 ${ISSM_DIR}/src/m/io/loadmodel.js \
    138                                 ${ISSM_DIR}/src/m/io/writetofile.js \
    139                                 ${ISSM_DIR}/src/m/io/writetofile.js \
     138                                ${ISSM_DIR}/src/m/io/*.js \
    140139                                ${ISSM_DIR}/src/m/materials/paterson.js \
    141140                                ${ISSM_DIR}/src/m/mesh/triangle.js \
    … …  
    145144                                ${ISSM_DIR}/src/m/parameterization/setflowequation.js \
    146145                                ${ISSM_DIR}/src/m/parameterization/setmask.js \
    147                                 ${ISSM_DIR}/src/m/plot/slider.js \
    148                                 ${ISSM_DIR}/src/m/plot/applyoptions.js \
    149                                 ${ISSM_DIR}/src/m/plot/checkplotoptions.js \
    150                                 ${ISSM_DIR}/src/m/plot/plot_manager.js \
    151                                 ${ISSM_DIR}/src/m/plot/plot_mesh.js \
    152                                 ${ISSM_DIR}/src/m/plot/plot_overlay.js \
    153                                 ${ISSM_DIR}/src/m/plot/plot_unit.js \
    154                                 ${ISSM_DIR}/src/m/plot/plot_quiver.js \
    155                                 ${ISSM_DIR}/src/m/plot/plotdoc.js \
    156                                 ${ISSM_DIR}/src/m/plot/plotmodel.js \
    157                                 ${ISSM_DIR}/src/m/plot/processdata.js \
    158                                 ${ISSM_DIR}/src/m/plot/processmesh.js \
    159                                 ${ISSM_DIR}/src/m/plot/gauge.js \
    160                                 ${ISSM_DIR}/src/m/plot/webgl.js \
     146                                ${ISSM_DIR}/src/m/plot/*.js \
    161147                                ${ISSM_DIR}/src/m/print/sprintf.js \
    162148                                ${ISSM_DIR}/src/m/solve/loadresultsfrombuffer.js \

  • issm/trunk-jpl/src/m/array/arrayoperations.js

    r21911 r22098  
    9393        return or;
    9494
     95} //}}}
     96function ArrayTranspose(array) { //{{{
     97        return Object.keys(array[0]).map(function(column) {
     98                return array.map(function(row) {
     99                        return row[column];
     100                });
     101        });
     102} //}}}
     103function ArrayConcat(a, b) { //{{{
     104        // Make sure that both typed arrays are of the same type
     105        if(Object.prototype.toString.call(a) !== Object.prototype.toString.call(b)) {
     106                throw 'The types of the two arguments passed for parameters a and b do not match.';
     107        }
     108        var c;
     109        if(a instanceof Float64Array) {
     110                c = new a.constructor(a.length + b.length);
     111                c.set(a);
     112                c.set(b, a.length);
     113        }
     114        else {
     115                c = clone(a).concat(clone(b));
     116        }       
     117        return c;
    95118} //}}}
    96119function ListToMatrix(list, elementsPerSubArray) { //{{{
    … …  
    181204        for(var i=0;i<array.length;i++){
    182205            for(var j=0;j<array[0].length;j++){
    183                 if (isNaN(array[i][j])) return 1;
     206                if (Number.isNaN(array[i][j])) return 1;
    184207            }
    185208        }
    … …  
    187210    else{
    188211        for(var i=0;i<array.length;i++){
    189             if (isNaN(array[i])) return 1;
     212            if (Number.isNaN(array[i])) return 1;
    190213        }
    191214    }
    … …  
    207230function ArrayAnyEqual(array,value) { //{{{
    208231       
    209         if(!isNaN(value)){
     232        if(!Number.isNaN(value)){
    210233                for(var i=0;i<array.length;i++){
    211234                        if (array[i]==value)return 1;
    … …  
    214237        else{
    215238                for(var i=0;i<array.length;i++){
    216                         if (isNaN(array[i]))return 1;
     239                        if (Number.isNaN(array[i]))return 1;
    217240                }
    218241        }
    … …  
    267290} //}}}
    268291function NewArrayFill(size,value) { //{{{
    269        
    270292        var array = new Array(size);
    271293       
    … …  
    276298        return array;
    277299} //}}}
     300function NewArrayFill2D(rows,cols,value) { //{{{
     301        var arr=new Array(rows);
     302
     303        for (var i=0;i<rows;i++) {
     304                arr[i] = NewArrayFill(cols,value);
     305        }
     306
     307        return arr;
     308} //}}}
    278309function NewArrayFillIncrement(size,start,increment) { //{{{
    279 
    280310        var array=new Array(size);
    281311
    … …  
    287317} //}}}
    288318function ArrayFind(array,value) { //{{{
    289        
    290         //find number of indices
    291319        var count=0;
    292         for (var i=0;i<array.length;i++)if(array[i]==value)count++;
    293 
    294         //allocate:
    295         var indices= NewArrayFill(count,0);
    296 
    297         //fill in:
    298         count=0;
    299         for (var i=0;i<array.length;i++){
    300                 if(array[i]==value){
    301                         indices[count]=i;
     320        var indices=[];
     321
     322        for (var i=0;i<array.length;i++){
     323                if(array[i]===value){
     324                        indices.push(count);
     325                }
     326                count++;
     327        }
     328        return indices;
     329} //}}}
     330function ArrayFind2D(array,value) { //{{{
     331        var count=0;
     332        var indices=[];
     333
     334        for (var i=0;i<array.length;i++){
     335                for (var j=0;j<array[i].length;j++){
     336                        if(array[i][j]===value){
     337                                indices.push(count);
     338                        }
    302339                        count++;
    303340                }
    … …  
    306343} //}}}
    307344function ArrayFindNot(array,value) { //{{{
    308        
    309         //find number of indices
    310345        var count=0;
    311         for (var i=0;i<array.length;i++)if(array[i]!=value)count++;
    312 
    313         //allocate:
    314         var indices= NewArrayFill(count,0);
    315 
    316         //fill in:
    317         count=0;
    318         for (var i=0;i<array.length;i++){
    319                 if(array[i]!=value){
    320                         indices[count]=i;
     346        var indices=[];
     347
     348        for (var i=0;i<array.length;i++){
     349                if(array[i]!==value){
     350                        indices.push(count);
     351                }
     352                count++;
     353        }
     354        return indices;
     355} //}}}
     356function ArrayFindNot2D(array,value) { //{{{
     357        var count=0;
     358        var indices=[];
     359
     360        for (var i=0;i<array.length;i++){
     361                for (var j=0;j<array[i].length;j++){
     362                        if(array[i][j]!==value){
     363                                indices.push(count);
     364                        }
    321365                        count++;
    322366                }
    … …  
    325369} //}}}
    326370function Create2DArray(rows,cols) { //{{{
    327         var arr = [];
     371        var arr=new Array(rows);
    328372
    329373        for (var i=0;i<rows;i++) {

  • issm/trunk-jpl/src/m/classes/inversion.js

    r21065 r22098  
    8181        this.thickness_obs=project3d(md,'vector',this.thickness_obs,'type','node');
    8282
    83         if (numel(this.cost_functions_coefficients)>1) {
     83        if (this.cost_functions_coefficients.length>1) {
    8484            this.cost_functions_coefficients=project3d(md,'vector',this.cost_functions_coefficients,'type','node');
    8585        }                       
    86         if (numel(this.min_parameters)>1) {
     86        if (this.min_parameters.length>1) {
    8787            this.min_parameters=project3d(md,'vector',this.min_parameters,'type','node');
    8888        }
    89         if (numel(this.max_parameters)>1) {
     89        if (this.max_parameters.length>1) {
    9090            this.max_parameters=project3d(md,'vector',this.max_parameters,'type','node');
    9191        }

  • issm/trunk-jpl/src/m/classes/model.js

    r21945 r22098  
    3434                        console.log(sprintf("//19s: //-22s -- //s","thermal"         ,"[1x1 " + typeof(this.thermal) + "]","parameters for thermal solution"));
    3535                        console.log(sprintf("//19s: //-22s -- //s","steadystate"     ,"[1x1 " + typeof(this.steadystate) + "]","parameters for steadystate solution"));
    36                         console.log(sprintf("//19s: //-22s -- //s","trans"       ,"[1x1 " + typeof(this.trans) + "]","parameters for trans solution"));
     36                        console.log(sprintf("//19s: //-22s -- //s","trans"           ,"[1x1 " + typeof(this.trans) + "]","parameters for trans solution"));
    3737                        console.log(sprintf("//19s: //-22s -- //s","levelset"        ,"[1x1 " + typeof(this.levelset) + "]","parameters for moving boundaries (level-set method)"));
    3838                        console.log(sprintf("//19s: //-22s -- //s","calving"         ,"[1x1 " + typeof(this.calving) + "]","parameters for calving"));
    3939                        console.log(sprintf("//19s: //-22s -- //s","gia"             ,"[1x1 " + typeof(this.gia) + "]","parameters for gia solution"));
     40                        console.log(sprintf("//19s: //-22s -- //s","love"            ,"[1x1 " + typeof(this.love) + "]","parameters for love solution"));
    4041                        console.log(sprintf("//19s: //-22s -- //s","esa"             ,"[1x1 " + typeof(this.esa) + "]","parameters for elastic adjustment solution"));
    4142                        console.log(sprintf("//19s: //-22s -- //s","autodiff"        ,"[1x1 " + typeof(this.autodiff) + "]","automatic differentiation parameters"));
    … …  
    5556                        this.constants        = new constants();
    5657                        this.geometry         = new geometry();
    57                         this.initialization   = new initialization();
    5858                        this.smb              = new SMBforcing();
    5959                        this.basalforcings    = new basalforcings();
     60                        this.materials        = new matice();
     61                        this.damage           = new damage();
    6062                        this.friction         = new friction();
     63                        this.flowequation     = new flowequation();
     64                        this.timestepping     = new timestepping();
     65                        this.initialization   = new initialization();
    6166                        this.rifts            = new rifts();
    6267                        this.slr              = new slr();
    63                         this.timestepping     = new timestepping();
    64                         this.groundingline    = new groundingline();
    65                         this.materials        = new matice();
    66                         this.damage           = new damage();
    67                         this.flowequation     = new flowequation();
     68
    6869                        this.debug            = new debug();
    6970                        this.verbose          = new verbose();
    … …  
    7172                        this.toolkits         = new toolkits();
    7273                        this.cluster          = new local();
     74
    7375                        this.balancethickness = new balancethickness();
    7476                        this.stressbalance    = new stressbalance();
     77                        this.groundingline    = new groundingline();
    7578                        this.hydrology        = new hydrologyshreve();
    7679                        this.masstransport    = new masstransport();
    … …  
    8083                        this.levelset         = new levelset();
    8184                        this.calving          = new calving();
     85                        this.love             = new fourierlove();
    8286                        this.gia              = new giaivins();
     87
    8388                        this.esa              = new esa();
    8489                        this.autodiff         = new autodiff();
    … …  
    8691                        this.qmu              = new qmu();
    8792                        this.amr              = new amr();
    88                         this.radaroverlay     = new radaroverlay();
     93
    8994                        this.results          = {};
    9095                        this.outputdefinition = new outputdefinition();
     96                        this.radaroverlay     = new radaroverlay();
    9197                        this.miscellaneous    = new miscellaneous();
    9298                        this.priv             = new priv();
    … …  
    113119
    114120                } //}}}
    115         this.extrude = function(md) { //{{{
    116         //EXTRUDE - vertically extrude a 2d mesh
    117         //
    118         //   vertically extrude a 2d mesh and create corresponding 3d mesh.
    119         //   The vertical distribution can:
    120         //    - follow a polynomial law
    121         //    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
    122         //    - be discribed by a list of coefficients (between 0 and 1)
    123         //   
    124         //
    125         //   Usage:
    126         //      md=extrude(md,numlayers,extrusionexponent);
    127         //      md=extrude(md,numlayers,lowerexponent,upperexponent);
    128         //      md=extrude(md,listofcoefficients);
    129         //
    130         //   Example:
    131         //      md=extrude(md,15,1.3);
    132         //      md=extrude(md,15,1.3,1.2);
    133         //      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
    134         //
    135         //   See also: MODELEXTRACT, COLLAPSE
    136 
    137         //some checks on list of arguments
    138             var argc = arguments.length;
    139 
    140             if ((argc > 4) | (argc < 2)) {
    141                 throw "extrude error message";
    142             }
    143 
    144             //Extrude the mesh
    145             if (argc==2) { //list of coefficients
    146                 clist=arguments[0];
    147 
    148                 if (clist.some(function(x) { return x < 0; }) | clist.some(function(x) { return x > 1; }))
    149                     throw 'extrusioncoefficients must be between 0 and 1';
    150                 clistcopy.push.apply(clistcopy, [0, 1]);
    151 
    152                 extrusionlist=clistcopy.reverse().filter(function(e, i, arr) { //only keep unique elements
    153                     return arr.indexOf(e, i+1) === -1;
    154                 }).sort();
    155 
    156                 numlayers=extrusionlist.length;
    157             } else if (argc==3) { //one polynomial law
    158                 if (arguments[1]<=0) {
    159                     //help extrude;
    160                     throw 'extrusionexponent must be >=0';
    161                 }
    162                 numlayers=arguments[0];
    163 
    164                 extrusionlist = [];
    165 
    166                 for (var i = 0; i < numlayers; i++) {
    167                     extrusionlist.push(Math.pow(i/(numlayers-1), arguments[1]));
    168                 };
    169 
    170             } else if (argc==4) { //two polynomial laws
    171                 numlayers=arguments[0];
    172                 lowerexp=arguments[1];
    173                 upperexp=arguments[2];
    174 
    175                 if (arguments[1]<=0 | args[2]<=0) {
    176                     throw 'lower and upper extrusionexponents must be >=0';
    177                 }
    178 
    179                 lowerextrusionlist = [];
    180                 for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
    181                     lowerextrusionlist.push(Math.pow(i, lowerexp)/2)
    182                 };
    183 
    184                 upperextrusionlist = [];
    185                 for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
    186                     upperextrusionlist.push(Math.pow(i, upperexp)/2)
    187                 };
    188 
    189                 //lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2;
    190                 //upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2;
    191 
    192                 var temp = lowerextrusionlist;
    193                 temp.push.apply(temp, [1]);
    194 
    195                 var temp2 = upperextrusionlist.map(function(x) { return -x; });
    196                 temp.push.apply(temp, temp2);
    197 
    198                 clistcopy.push.apply(clistcopy, [0, 1]);
    199                 //extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist]));
    200                 extrusionlist=temp.sort();
    201             }
    202 
    203             if (numlayers<2) {
    204                 console.error('number of layers should be at least 2');
    205             }
    206 
    207             if (md.mesh.domaintype() === '3D') {
    208                 console.error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
    209             }
    210 
    211             //Initialize with the 2d mesh
    212             var mesh2d = md.mesh;
    213             md.mesh=new mesh3dprisms();
    214 
    215             md.mesh.x                           = mesh2d.x;
    216             md.mesh.y                           = mesh2d.y;
    217             md.mesh.elements                    = mesh2d.elements;
    218             md.mesh.numberofelements            = mesh2d.numberofelements;
    219             md.mesh.numberofvertices            = mesh2d.numberofvertices;
    220 
    221             md.mesh.lat                         = mesh2d.lat;
    222             md.mesh.long                        = mesh2d.long;
    223             md.mesh.epsg                        = mesh2d.epsg;
    224 
    225             md.mesh.vertexonboundary            = mesh2d.vertexonboundary;
    226             md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity;
    227             md.mesh.elementconnectivity         = mesh2d.elementconnectivity;
    228             md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
    229 
    230             md.mesh.extractedvertices           = mesh2d.extractedvertices;
    231             md.mesh.extractedelements           = mesh2d.extractedelements;
    232 
    233             x3d=[];
    234             y3d=[];
    235             z3d=[];  //the lower node is on the bed
    236             thickness3d=md.geometry.thickness; //thickness and bed for these nodes
    237                 bed3d=md.geometry.base;
    238 
    239             //Create the new layers
    240             for (var i = 1; i <= numlayers; i++) {
    241                 x3d=[x3d, md.mesh.x];
    242                 y3d=[y3d, md.mesh.y];
    243                 //nodes are distributed between bed and surface accordingly to the given exponent
    244                 z3d=[z3d, bed3d+thickness3d*extrusionlist(i)];
    245             };
    246 
    247             number_nodes3d=x3d.length; //number of 3d nodes for the non extruded part of the mesh
    248 
    249             //Extrude elements
    250             elements3d=[];
    251             for (var i = 0; i < numlayers; i++) {
    252                 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
    253             };
    254 
    255             number_el3d=elements3d.length; //number of 3d nodes for the non extruded part of the mesh
    256 
    257             function fillArray(value, len) {
    258                 if (len == 0) return [];
    259                 var a = [value];
    260                 while (a.length * 2 <= len) a = a.concat(a);
    261                 if (a.length < len) a = a.concat(a.slice(0, len - a.length));
    262                 return a;
    263             }
    264 
    265             //Keep a trace of lower and upper nodes
    266             //lowervertex=NaN*ones(number_nodes3d,1);
    267             //uppervertex=NaN*ones(number_nodes3d,1);
    268 
    269             lowervertex = fillArray(NaN, number_nodes3d);
    270             uppervertex = fillArray(NaN, number_nodes3d);
    271 
    272             //lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices;
    273             //uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d;
    274 
    275             for (var i = md.mesh.numberofvertices+1, k = 1; i < lowervertex.length && k <= (numlayers-1)*md.mesh.numberofvertices; i++, k++) {
    276                 lowervertex[i] = k;
    277             };
    278 
    279             for (var i = 1, k = md.mesh.numberofvertices+1; i <= (numlayers-1)*md.mesh.numberofvertices && k <= number_nodes3d; i++, k++) {
    280                 uppervertex[i] = k;
    281             };
    282 
    283             md.mesh.lowervertex=lowervertex;
    284             md.mesh.uppervertex=uppervertex;
    285 
    286             //same for lower and upper elements
    287             lowerelements = fillArray(NaN, number_el3d);
    288             upperelements = fillArray(NaN, number_el3d);
    289             //lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements;
    290             //upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements;
    291 
    292             for (var i = md.mesh.numberofvertices+1, k = 1; i < lowerelements.length && k <= (numlayers-2)*md.mesh.numberofelements; i++, k++) {
    293                 lowerelements[i] = k;
    294             };
    295 
    296             for (var i = 1, k = md.mesh.numberofelements + 1; i <= (numlayers-2)*md.mesh.numberofelements && k <= (numlayers-1)*md.mesh.numberofelements; i++, k++) {
    297                 upperelements[i] = k;
    298             };
    299 
    300             md.mesh.lowerelements=lowerelements;
    301             md.mesh.upperelements=upperelements;
    302 
    303             //Save old mesh
    304             md.mesh.x2d=md.mesh.x;
    305             md.mesh.y2d=md.mesh.y;
    306             md.mesh.elements2d=md.mesh.elements;
    307             md.mesh.numberofelements2d=md.mesh.numberofelements;
    308             md.mesh.numberofvertices2d=md.mesh.numberofvertices;
    309 
    310             //Build global 3d mesh
    311             md.mesh.elements=elements3d;
    312             md.mesh.x=x3d;
    313             md.mesh.y=y3d;
    314             md.mesh.z=z3d;
    315             md.mesh.numberofelements=number_el3d;
    316             md.mesh.numberofvertices=number_nodes3d;
    317             md.mesh.numberoflayers=numlayers;
    318 
    319             //Ok, now deal with the other fields from the 2d mesh:
    320 
    321             //bedinfo and surface info
    322             md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
    323             md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
    324             md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
    325 
    326             //lat long
    327             md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
    328             md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
    329 
    330                 //md.geometry=md.geometry.extrude(md);
    331                 //md.friction=md.friction.extrude(md);
    332                 //md.inversion=md.inversion.extrude(md);
    333                 //md.smb=md.smb.extrude(md);
    334                 //md.initialization=md.initialization.extrude(md);
    335 
    336                 //md.flowequation=md.flowequation.extrude(md);
    337                 //md.stressbalance=md.stressbalance.extrude(md);
    338                 //md.thermal=md.thermal.extrude(md);
    339                 //md.masstransport=md.masstransport.extrude(md);
    340                 //md.levelset=md.levelset.extrude(md);
    341                 //md.calving=md.calving.extrude(md);
    342                 //md.hydrology = md.hydrology.extrude(md);
    343 
    344                 ////connectivity
    345                 ////if ~isnan(md.mesh.elementconnectivity)
    346 
    347                 //if (md.mesh.elementconnectivity.every(function(e, i, arr) { return e !== null && typeof e !== 'undefined'; })) {
    348                     //var temparr = [];
    349                     //for (var i = 0; i < numlayers; i++) {
    350                         //temparr.push(md.mesh.elementconnectivity);
    351                     //};
    352 
    353                     ////md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
    354 
    355                     //md.mesh.elementconnectivity = temparr;
    356 
    357                     ////md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
    358 
    359                     //for (var i = 0; i < md.mesh.elementconnectivity.length; i++) {
    360                         //if (md.mesh.elementconnectivity[i] == 0) {
    361                             //md.mesh.elementconnectivity[i] = NaN;
    362                         //}
    363                     //};
    364 
    365                     //for (var i = 2; i < numlayers; i++) {
    366                         //for (var j = (i-1)*md.mesh.nuberofelements2d+1; i <= (i)*md.mesh.numberofelements2d; i++) {
    367                             //md.mesh.elementconnectivity[i] += md.mesh.numberofelements2d;
    368                         //};
    369                         ////md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
    370                             ////=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
    371                     //};
    372 
    373                     //md.mesh.elementconnectivity.map(function(x) { return (isNaN(x)) ? 0 : x; });
    374                 //}
    375 
    376                 //md.materials=md.materials.extrude(md);
    377                 //md.damage=md.damage.extrude(md);
    378                 //md.mask=md.mask.extrude(md);
    379                 //md.qmu=md.qmu.extrude(md);
    380                 //md.basalforcings=md.basalforcings.extrude(md);
    381 
    382                 ////increase connectivity if less than 25:
    383                 //if (md.mesh.average_vertex_connectivity<=25)
    384                     //md.mesh.average_vertex_connectivity=100;
    385         } //}}}
    386         this.collapse = function(md) { //{{{
    387             /*
    388              *COLLAPSE - collapses a 3d mesh into a 2d mesh
    389              *
    390              *   This routine collapses a 3d model into a 2d model
    391              *   and collapses all the fileds of the 3d model by
    392              *   taking their depth-averaged values
    393              *
    394              *   Usage:
    395              *     md=collapse(md)
    396              *
    397              *   See also: EXTRUDE, MODELEXTRACT
    398              */
    399 
    400             // Check that the model is really a 3d model
    401             if (md.mesh.elementtype() !== 'Penta') {
    402                 console.error('collapse error message: only 3d mesh can be collapsed')
    403             }
    404 
    405             // Start with changing all the fields from the 3d mesh
    406 
    407             // dealing with the friction law
    408             // drag is limited to nodes that are on the bedrock.
    409             if (md.friction.classname() === 'friction') {
    410                 md.friction.coefficient=project2d(md,md.friction.coefficient,1);
    411                 md.friction.p=project2d(md,md.friction.p,1);
    412                 md.friction.q=project2d(md,md.friction.q,1);
    413             } else if (md.friction.classname() === 'frictionhydro') {
    414                 md.friction.q=project2d(md,md.friction.q,1);
    415                 md.friction.C=project2d(md,md.friction.C,1);
    416                 md.friction.As=project2d(md,md.friction.As,1);
    417                 md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
    418             } else if (md.friction.classname() === 'frictionwaterlayer') {
    419                 md.friction.coefficient=project2d(md,md.friction.coefficient,1);
    420                 md.friction.p=project2d(md,md.friction.p,1);
    421                 md.friction.q=project2d(md,md.friction.q,1);
    422                 md.friction.water_layer=project2d(md,md.friction.water_layer,1);
    423             } else if (md.friction.classname() === 'frictionweertman') {
    424                 md.friction.C=project2d(md,md.friction.C,1);
    425                 md.friction.m=project2d(md,md.friction.m,1);
    426             } else if (md.friction.classname() === 'frictionweertmantemp') {
    427                 md.friction.C=project2d(md,md.friction.C,1);
    428                 md.friction.m=project2d(md,md.friction.m,1);
    429             } else {
    430                 disp('friction type not supported');
    431             }
    432 
    433             // observations
    434             if (!isNaN(md.inversion.vx_obs))
    435                 md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
    436 
    437             if (!isNaN(md.inversion.vy_obs))
    438                 md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
    439 
    440             if (!isNaN(md.inversion.vel_obs))
    441                 md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
    442 
    443             if (!isNaN(md.inversion.cost_functions_coefficients))
    444                 md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
    445 
    446             if (numel(md.inversion.min_parameters)>1)
    447                 md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
    448 
    449             if (numel(md.inversion.max_parameters)>1)
    450                 md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
    451 
    452             if (md.smb.classname() === 'SMBforcing' && !isNaN(md.smb.mass_balance)) {
    453                 md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
    454             } else if (md.smb.classname() === 'SMBhenning' && !isNaN(md.smb.smbref)) {
    455                 md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
    456             }
    457 
    458             // results
    459             if (!isNaN(md.initialization.vx))
    460                 md.initialization.vx=DepthAverage(md,md.initialization.vx);
    461             if (!isNaN(md.initialization.vy))
    462                 md.initialization.vy=DepthAverage(md,md.initialization.vy);
    463             if (!isNaN(md.initialization.vz))
    464                 md.initialization.vz=DepthAverage(md,md.initialization.vz);
    465             if (!isNaN(md.initialization.vel))
    466                 md.initialization.vel=DepthAverage(md,md.initialization.vel);
    467             if (!isNaN(md.initialization.temperature))
    468                 md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
    469             if (!isNaN(md.initialization.pressure))
    470                 md.initialization.pressure=project2d(md,md.initialization.pressure,1);
    471             if (!isNaN(md.initialization.sediment_head))
    472                 md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
    473             if (!isNaN(md.initialization.epl_head))
    474                 md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
    475             if (!isNaN(md.initialization.epl_thickness))
    476                 md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
    477 
    478             // giaivins
    479             if (!isNaN(md.gia.mantle_viscosity))
    480                 md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1);
    481             if (!isNaN(md.gia.lithosphere_thickness))
    482                 md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1);
    483 
    484             // elementstype
    485             if (!isNaN(md.flowequation.element_equation)) {
    486                 md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
    487                 md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
    488                 md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
    489                 md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
    490                 md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
    491             }
    492 
    493             // boundary conditions
    494             md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
    495             md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
    496             md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
    497             md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
    498             md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
    499             md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
    500             if (!isNaN(md.damage.spcdamage))
    501                 md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
    502             md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
    503 
    504             // Hydrologydc variables
    505             if (md.hydrology.classname() === 'hydrologydc') {
    506                 md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
    507                 md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
    508                 md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
    509                 md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
    510 
    511                 if(md.hydrology.isefficientlayer==1)
    512                     md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
    513             }
    514            
    515             // materials
    516             md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
    517             md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
    518            
    519             // damage:
    520             if (md.damage.isdamage)
    521                 md.damage.D=DepthAverage(md,md.damage.D);
    522 
    523             // special for thermal modeling:
    524             if (!isNaN(md.basalforcings.groundedice_melting_rate))
    525                 md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1);
    526 
    527             if (!isNaN(md.basalforcings.floatingice_melting_rate))
    528                 md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1);
    529 
    530             md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); // bedrock only gets geothermal flux
    531 
    532             // update of connectivity matrix
    533             md.mesh.average_vertex_connectivity=25;
    534 
    535             // Collapse the mesh
    536             var nodes2d=md.mesh.numberofvertices2d;
    537             var elements2d=md.mesh.numberofelements2d;
    538 
    539             // parameters
    540             md.geometry.surface=project2d(md,md.geometry.surface,1);
    541             md.geometry.thickness=project2d(md,md.geometry.thickness,1);
    542             md.geometry.base=project2d(md,md.geometry.base,1);
    543             if (!isNaN(md.geometry.bed))
    544                 md.geometry.bed=project2d(md,md.geometry.bed,1);
    545 
    546             if (!isNaN(md.mask.groundedice_levelset))
    547                 md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1);
    548 
    549             if (!isNaN(md.mask.ice_levelset))
    550                 md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
    551 
    552             // lat long
    553             if (numel(md.mesh.lat) === md.mesh.numberofvertices)
    554                 md.mesh.lat=project2d(md,md.mesh.lat,1);
    555             if (numel(md.mesh.long) === md.mesh.numberofvertices)
    556                 md.mesh.long=project2d(md,md.mesh.long,1);
    557 
    558             // Initialize with the 2d mesh
    559             var mesh = new mesh2d();
    560             mesh.x=md.mesh.x2d;
    561             mesh.y=md.mesh.y2d;
    562             mesh.numberofvertices=md.mesh.numberofvertices2d;
    563             mesh.numberofelements=md.mesh.numberofelements2d;
    564             mesh.elements=md.mesh.elements2d;
    565 
    566             if (!isNaN(md.mesh.vertexonboundary))
    567                 mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
    568             if (!isNaN(md.mesh.elementconnectivity))
    569                 mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
    570 
    571             md.mesh=mesh;
    572             md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
    573             md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
    574             md.mesh.segments=contourenvelope(md.mesh);
    575 
    576             return md;
    577         } /*}}}*/
     121                this.extrude = function(md) { //{{{
     122                        //EXTRUDE - vertically extrude a 2d mesh
     123                        //
     124                        //   vertically extrude a 2d mesh and create corresponding 3d mesh.
     125                        //   The vertical distribution can:
     126                        //    - follow a polynomial law
     127                        //    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
     128                        //    - be described by a list of coefficients (between 0 and 1)
     129                        //   
     130                        //
     131                        //   Usage:
     132                        //      md=extrude(md,numlayers,extrusionexponent);
     133                        //      md=extrude(md,numlayers,lowerexponent,upperexponent);
     134                        //      md=extrude(md,listofcoefficients);
     135                        //
     136                        //   Example:
     137                        //      md=extrude(md,15,1.3);
     138                        //      md=extrude(md,15,1.3,1.2);
     139                        //      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
     140                        //
     141                        //   See also: MODELEXTRACT, COLLAPSE
     142
     143                        //some checks on list of arguments
     144                        var argc = arguments.length;
     145                        var extrusionlist;
     146                        var numlayers;
     147                       
     148                        if ((argc > 4) | (argc < 2)) {
     149                                throw "extrude error message";
     150                        }
     151
     152                        //Extrude the mesh
     153                        if (argc==2) { //list of coefficients
     154                                clist=arguments[0];
     155                                if (ArrayAnyBelowStrict(clist,0) || ArrayAnyAboveStrict(clist,1)) {
     156                                        throw 'extrusioncoefficients must be between 0 and 1';
     157                                }
     158                                extrusionlist=ArraySort(ArrayUnique(clist.push(0,1)));
     159                                numlayers=extrusionlist.length;
     160                        } else if (argc==3) { //one polynomial law
     161                                if (arguments[1]<=0) {
     162                                        throw 'extrusionexponent must be >=0';
     163                                }
     164                                numlayers=arguments[1];
     165                                extrusionlist = [];
     166                                for (var i = 0; i < numlayers; i++) {
     167                                        extrusionlist.push(Math.pow(i/(numlayers-1), arguments[2]));
     168                                }
     169                        } else if (argc==4) { //two polynomial laws
     170                                numlayers=arguments[1];
     171                                var lowerexp=arguments[2];
     172                                var upperexp=arguments[3];
     173
     174                                if (arguments[2]<=0 || arguments[3]<=0) {
     175                                        throw 'lower and upper extrusionexponents must be >=0';
     176                                }
     177
     178                                var lowerextrusionlist = [];
     179                                for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
     180                                        lowerextrusionlist.push(Math.pow(i, lowerexp)/2)
     181                                }
     182                                var upperextrusionlist = [];
     183                                for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
     184                                        upperextrusionlist.push(1-Math.pow(i, upperexp)/2)
     185                                }
     186                                extrusionlist=ArrayConcat(lowerextrusionlist,upperextrusionlist);
     187                                extrusionlist=ArraySort(ArrayUnique(extrusionlist.push(1)));
     188                        }
     189
     190                        if (numlayers<2) {
     191                                console.error('number of layers should be at least 2');
     192                        }
     193                        if (md.mesh.domaintype() === '3D') {
     194                                console.error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
     195                        }
     196
     197                        //Initialize with the 2d mesh
     198                        var mesh2d = md.mesh;
     199                        md.mesh=new mesh3dprisms();
     200                        md.mesh.x                           = mesh2d.x;
     201                        md.mesh.y                           = mesh2d.y;
     202                        md.mesh.elements                    = mesh2d.elements;
     203                        md.mesh.numberofelements            = mesh2d.numberofelements;
     204                        md.mesh.numberofvertices            = mesh2d.numberofvertices;
     205
     206                        md.mesh.lat                         = mesh2d.lat;
     207                        md.mesh.long                        = mesh2d.long;
     208                        md.mesh.epsg                        = mesh2d.epsg;
     209
     210                        md.mesh.vertexonboundary            = mesh2d.vertexonboundary;
     211                        md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity;
     212                        md.mesh.elementconnectivity         = mesh2d.elementconnectivity;
     213                        md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
     214
     215                        md.mesh.extractedvertices           = mesh2d.extractedvertices;
     216                        md.mesh.extractedelements           = mesh2d.extractedelements;
     217
     218                        var x3d=new Float64Array();
     219                        var y3d=new Float64Array();
     220                        var z3d=new Float64Array();  //the lower node is on the bed
     221                        var thickness3d=md.geometry.thickness; //thickness and bed for these nodes
     222                        var bed3d=md.geometry.base;
     223
     224                        //Create the new layers
     225                        //Dimensions of x/y/z3d: md.mesh.numberofvertices * numlayers, 1
     226                        for (var i = 1; i <= numlayers; i++) {
     227                                x3d=ArrayConcat(x3d, md.mesh.x);
     228                                y3d=ArrayConcat(y3d, md.mesh.y);
     229                                //nodes are distributed between bed and surface accordingly to the given exponent
     230                                z3d=ArrayConcat(z3d, FloatFix(thickness3d.map(function(value,index) { return bed3d[index] + value * extrusionlist[i-1]; }), thickness3d.length));
     231                        }
     232                        var number_nodes3d=x3d.length; //number of 3d nodes for the non extruded part of the mesh
     233
     234                        //Extrude elements
     235                        //Create the elements of the 3d mesh for the non extruded part
     236                        //Dimensions of elements3d: md.mesh.numberofelements * (numlayers - 1), 6
     237                        var elements3d=[];
     238                        var elements_prisms=[];
     239                        for (var i = 1; i < numlayers; i++) {
     240                                for (var j = 0; j < md.mesh.numberofelements; j++) {
     241                                        elements_prisms = [];
     242                                        for (var k = 0, offset = (i - 1) * md.mesh.numberofvertices; k < 3; k++) {
     243                                                elements_prisms.push(md.mesh.elements[j][k]+offset);
     244                                        }
     245                                        for (var k = 0, offset = i * md.mesh.numberofvertices; k < 3; k++) {
     246                                                elements_prisms.push(md.mesh.elements[j][k]+offset);
     247                                        }
     248                                        elements3d.push(elements_prisms);
     249                                }                               
     250                        }
     251                        number_el3d=elements3d.length; //number of 3d nodes for the non extruded part of the mesh
     252
     253                        //Keep a trace of lower and upper nodes
     254                        var lowervertex = NewArrayFill(number_nodes3d, NaN);
     255                        var uppervertex = NewArrayFill(number_nodes3d, NaN);
     256                        //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofvertices+1)
     257                        for (var i = md.mesh.numberofvertices, k = 1; i < lowervertex.length; i++, k++) {
     258                                lowervertex[i] = k;
     259                        };
     260                        //Set last layer to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1)
     261                        for (var i = 0, k = md.mesh.numberofvertices+1; i < (numlayers-1)*md.mesh.numberofvertices; i++, k++) {
     262                                uppervertex[i] = k;
     263                        };
     264                        md.mesh.lowervertex=lowervertex;
     265                        md.mesh.uppervertex=uppervertex;
     266
     267                        //same for lower and upper elements
     268                        var lowerelements = NewArrayFill(number_el3d, NaN);
     269                        var upperelements = NewArrayFill(number_el3d, NaN);
     270                        //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofelements+1)
     271                        for (var i = md.mesh.numberofelements, k = 1; i < lowerelements.length; i++, k++) {
     272                                lowerelements[i] = k;
     273                        };
     274                        //Set last 2 layers to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1)
     275                        for (var i = 0, k = md.mesh.numberofelements + 1; i < (numlayers-2)*md.mesh.numberofelements; i++, k++) {
     276                                upperelements[i] = k;
     277                        };
     278                        md.mesh.lowerelements=lowerelements;
     279                        md.mesh.upperelements=upperelements;
     280
     281                        //Save old mesh
     282                        md.mesh.x2d=md.mesh.x;
     283                        md.mesh.y2d=md.mesh.y;
     284                        md.mesh.elements2d=md.mesh.elements;
     285                        md.mesh.numberofelements2d=md.mesh.numberofelements;
     286                        md.mesh.numberofvertices2d=md.mesh.numberofvertices;
     287
     288                        //Build global 3d mesh
     289                        md.mesh.elements=elements3d;
     290                        md.mesh.x=x3d;
     291                        md.mesh.y=y3d;
     292                        md.mesh.z=z3d;
     293                        md.mesh.numberofelements=number_el3d;
     294                        md.mesh.numberofvertices=number_nodes3d;
     295                        md.mesh.numberoflayers=numlayers;
     296
     297                        //Ok, now deal with the other fields from the 2d mesh:
     298
     299                        //bedinfo and surface info
     300                        md.mesh.vertexonbase=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',1);
     301                        md.mesh.vertexonsurface=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
     302                        md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
     303
     304                        //lat long
     305                        md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
     306                        md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
     307
     308                        md.geometry=md.geometry.extrude(md);
     309                        md.friction=md.friction.extrude(md);
     310                        md.inversion=md.inversion.extrude(md);
     311                        md.smb=md.smb.extrude(md);
     312                        md.initialization=md.initialization.extrude(md);
     313
     314                        md.flowequation=md.flowequation.extrude(md);
     315                        md.stressbalance=md.stressbalance.extrude(md);
     316                        md.thermal=md.thermal.extrude(md);
     317                        md.masstransport=md.masstransport.extrude(md);
     318                        md.levelset=md.levelset.extrude(md);
     319                        md.calving=md.calving.extrude(md);
     320                        md.hydrology = md.hydrology.extrude(md);
     321
     322                        //connectivity
     323                        if (!ArrayAnyNaN(md.mesh.elementconnectivity)) {
     324                                //md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); //Replicate the matrix across numlayers-1
     325                                var temparr = [];
     326                                for (var i = 0; i < numlayers - 1; i++) {
     327                                        temparr = ArrayConcat(temparr, md.mesh.elementconnectivity);
     328                                }
     329                                md.mesh.elementconnectivity = temparr;
     330               
     331                                //md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
     332                                var indices = ArrayFind(md.mesh.elementconnectivity, 0);
     333                                for (var i = 0; i < indices.length; i++) {
     334                                        md.mesh.elementconnectivity[i] = NaN;
     335                                };
     336               
     337                                for (var i = 2; i < numlayers; i++) {
     338                                        //md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
     339                                        //=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
     340                                        for (var j = (i-1)*md.mesh.numberofelements2d; j <= (i)*md.mesh.numberofelements2d-1; j++) {
     341                                                for (var k = 0; k < 3; k++) {
     342                                                        md.mesh.elementconnectivity[j][k] += md.mesh.numberofelements2d;
     343                                                }
     344                                        }
     345                                }
     346                                md.mesh.elementconnectivity = md.mesh.elementconnectivity.map(function(value) { return (Number.isNaN(value)) ? 0 : value; });
     347                        }
     348
     349                        md.materials=md.materials.extrude(md);
     350                        md.damage=md.damage.extrude(md);
     351                        md.mask=md.mask.extrude(md);
     352                        md.qmu=md.qmu.extrude(md);
     353                        md.basalforcings=md.basalforcings.extrude(md);
     354
     355                        //increase connectivity if less than 25:
     356                        if (md.mesh.average_vertex_connectivity<=25) {
     357                                md.mesh.average_vertex_connectivity=100;
     358                        }
     359                } //}}}
     360                this.collapse = function(md) { //{{{
     361                        /*
     362                         *COLLAPSE - collapses a 3d mesh into a 2d mesh
     363                         *
     364                         *   This routine collapses a 3d model into a 2d model
     365                         *   and collapses all the fileds of the 3d model by
     366                         *   taking their depth-averaged values
     367                         *
     368                         *   Usage:
     369                         *       md=collapse(md)
     370                         *
     371                         *   See also: EXTRUDE, MODELEXTRACT
     372                         */
     373
     374                        // Check that the model is really a 3d model
     375                        if (md.mesh.elementtype() !== 'Penta') {
     376                                console.error('collapse error message: only 3d mesh can be collapsed')
     377                        }
     378
     379                        // Start with changing all the fields from the 3d mesh
     380
     381                        // dealing with the friction law
     382                        // drag is limited to nodes that are on the bedrock.
     383                        if (md.friction.classname() === 'friction') {
     384                                md.friction.coefficient=project2d(md,md.friction.coefficient,1);
     385                                md.friction.p=project2d(md,md.friction.p,1);
     386                                md.friction.q=project2d(md,md.friction.q,1);
     387                        } else if (md.friction.classname() === 'frictionhydro') {
     388                                md.friction.q=project2d(md,md.friction.q,1);
     389                                md.friction.C=project2d(md,md.friction.C,1);
     390                                md.friction.As=project2d(md,md.friction.As,1);
     391                                md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
     392                        } else if (md.friction.classname() === 'frictionwaterlayer') {
     393                                md.friction.coefficient=project2d(md,md.friction.coefficient,1);
     394                                md.friction.p=project2d(md,md.friction.p,1);
     395                                md.friction.q=project2d(md,md.friction.q,1);
     396                                md.friction.water_layer=project2d(md,md.friction.water_layer,1);
     397                        } else if (md.friction.classname() === 'frictionweertman') {
     398                                md.friction.C=project2d(md,md.friction.C,1);
     399                                md.friction.m=project2d(md,md.friction.m,1);
     400                        } else if (md.friction.classname() === 'frictionweertmantemp') {
     401                                md.friction.C=project2d(md,md.friction.C,1);
     402                                md.friction.m=project2d(md,md.friction.m,1);
     403                        } else {
     404                                disp('friction type not supported');
     405                        }
     406
     407                        // observations
     408                        if (!Number.isNaN(md.inversion.vx_obs))
     409                                md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
     410
     411                        if (!Number.isNaN(md.inversion.vy_obs))
     412                                md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
     413
     414                        if (!Number.isNaN(md.inversion.vel_obs))
     415                                md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
     416
     417                        if (!Number.isNaN(md.inversion.cost_functions_coefficients))
     418                                md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
     419
     420                        if (numel(md.inversion.min_parameters)>1)
     421                                md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
     422
     423                        if (numel(md.inversion.max_parameters)>1)
     424                                md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
     425
     426                        if (md.smb.classname() === 'SMBforcing' && !Number.isNaN(md.smb.mass_balance)) {
     427                                md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
     428                        } else if (md.smb.classname() === 'SMBhenning' && !Number.isNaN(md.smb.smbref)) {
     429                                md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
     430                        }
     431
     432                        // results
     433                        if (!Number.isNaN(md.initialization.vx))
     434                                md.initialization.vx=DepthAverage(md,md.initialization.vx);
     435                        if (!Number.isNaN(md.initialization.vy))
     436                                md.initialization.vy=DepthAverage(md,md.initialization.vy);
     437                        if (!Number.isNaN(md.initialization.vz))
     438                                md.initialization.vz=DepthAverage(md,md.initialization.vz);
     439                        if (!Number.isNaN(md.initialization.vel))
     440                                md.initialization.vel=DepthAverage(md,md.initialization.vel);
     441                        if (!Number.isNaN(md.initialization.temperature))
     442                                md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
     443                        if (!Number.isNaN(md.initialization.pressure))
     444                                md.initialization.pressure=project2d(md,md.initialization.pressure,1);
     445                        if (!Number.isNaN(md.initialization.sediment_head))
     446                                md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
     447                        if (!Number.isNaN(md.initialization.epl_head))
     448                                md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
     449                        if (!Number.isNaN(md.initialization.epl_thickness))
     450                                md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
     451
     452                        // giaivins
     453                        if (!Number.isNaN(md.gia.mantle_viscosity))
     454                                md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1);
     455                        if (!Number.isNaN(md.gia.lithosphere_thickness))
     456                                md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1);
     457
     458                        // elementstype
     459                        if (!Number.isNaN(md.flowequation.element_equation)) {
     460                                md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
     461                                md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
     462                                md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
     463                                md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
     464                                md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
     465                        }
     466
     467                        // boundary conditions
     468                        md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
     469                        md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
     470                        md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
     471                        md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
     472                        md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
     473                        md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
     474                        if (!Number.isNaN(md.damage.spcdamage))
     475                                md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
     476                        md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
     477
     478                        // Hydrologydc variables
     479                        if (md.hydrology.classname() === 'hydrologydc') {
     480                                md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
     481                                md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
     482                                md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
     483                                md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
     484
     485                                if(md.hydrology.isefficientlayer==1)
     486                                        md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
     487                        }
     488                       
     489                        // materials
     490                        md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
     491                        md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
     492                       
     493                        // damage:
     494                        if (md.damage.isdamage)
     495                                md.damage.D=DepthAverage(md,md.damage.D);
     496
     497                        // special for thermal modeling:
     498                        if (!Number.isNaN(md.basalforcings.groundedice_melting_rate))
     499                                md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1);
     500
     501                        if (!Number.isNaN(md.basalforcings.floatingice_melting_rate))
     502                                md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1);
     503
     504                        md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); // bedrock only gets geothermal flux
     505
     506                        // update of connectivity matrix
     507                        md.mesh.average_vertex_connectivity=25;
     508
     509                        // Collapse the mesh
     510                        var nodes2d=md.mesh.numberofvertices2d;
     511                        var elements2d=md.mesh.numberofelements2d;
     512
     513                        // parameters
     514                        md.geometry.surface=project2d(md,md.geometry.surface,1);
     515                        md.geometry.thickness=project2d(md,md.geometry.thickness,1);
     516                        md.geometry.base=project2d(md,md.geometry.base,1);
     517                        if (!Number.isNaN(md.geometry.bed))
     518                                md.geometry.bed=project2d(md,md.geometry.bed,1);
     519
     520                        if (!Number.isNaN(md.mask.groundedice_levelset))
     521                                md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1);
     522
     523                        if (!Number.isNaN(md.mask.ice_levelset))
     524                                md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
     525
     526                        // lat long
     527                        if (numel(md.mesh.lat) === md.mesh.numberofvertices)
     528                                md.mesh.lat=project2d(md,md.mesh.lat,1);
     529                        if (numel(md.mesh.long) === md.mesh.numberofvertices)
     530                                md.mesh.long=project2d(md,md.mesh.long,1);
     531
     532                        // Initialize with the 2d mesh
     533                        var mesh = new mesh2d();
     534                        mesh.x=md.mesh.x2d;
     535                        mesh.y=md.mesh.y2d;
     536                        mesh.numberofvertices=md.mesh.numberofvertices2d;
     537                        mesh.numberofelements=md.mesh.numberofelements2d;
     538                        mesh.elements=md.mesh.elements2d;
     539
     540                        if (!Number.isNaN(md.mesh.vertexonboundary))
     541                                mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
     542                        if (!Number.isNaN(md.mesh.elementconnectivity))
     543                                mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
     544
     545                        md.mesh=mesh;
     546                        md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
     547                        md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
     548                        md.mesh.segments=contourenvelope(md.mesh);
     549
     550                        return md;
     551                } /*}}}*/
    578552        //properties
    579553        // {{{
    … …  
    611585                this.levelset         = 0;
    612586                this.calving          = 0;
     587                this.love             = 0;
    613588                this.gia              = 0;
    614589                this.esa              = 0;

  • issm/trunk-jpl/src/m/classes/qmu.js

    r20838 r22098  
    9191        }// }}}
    9292    this.extrude = function(md) {//{{{
    93         this.partition=project3d(md,'vector',math.transpose(this.partition),'type','node');
     93        if (!isNaN(this.partition)) this.partition=project3d(md,'vector',ArrayTranspose(this.partition),'type','node');
    9494        return this;
    9595    }//}}}

  • issm/trunk-jpl/src/m/classes/thermal.js

    r21065 r22098  
    5555        this.spctemperature=project3d(md,'vector',this.spctemperature,'type','node','layer',md.mesh.numberoflayers,'padding',NaN);
    5656        if (md.initialization.temperature.length===md.mesh.numberofvertices) {
    57             this.spctemperature = ArrayFill2D(Create2DArray(md.mesh.numberofvertices, 1), NaN);
    58             var pos=ArrayFindNot2D(md.mesh.vertexonsurface, 0);
    59             pos.forEach(function(x) { // impose observed temperature on surface
    60                 this.spctemperature[x] = md.initialization.temperature[x];
    61             });
     57            this.spctemperature = NewArrayFill(md.mesh.numberofvertices, NaN);
     58            var pos=ArrayFindNot(md.mesh.vertexonsurface, 0);
     59                        // impose observed temperature on surface
     60                        for (var i=0,posIndex=0,count=0;i<md.initialization.temperature.length;i++){
     61                                if(count===pos[posIndex]){
     62                                        this.spctemperature[i] = md.initialization.temperature[i][0];
     63                                        posIndex++;
     64                                }
     65                                count++;
     66                        }
     67//            this.spctemperature = NewArrayFill2D(md.mesh.numberofvertices, 1, NaN);
     68//            var pos=ArrayFindNot2D(md.mesh.vertexonsurface, 0);
     69//                      // impose observed temperature on surface
     70//                      for (var i=0,posIndex=0,count=0;i<md.initialization.temperature.length;i++){
     71//                              for (var j=0;j<md.initialization.temperature[i].length;j++){
     72//                                      if(count===pos[posIndex]){
     73//                                              this.spctemperature[i][j] = md.initialization.temperature[i][j];
     74//                                              posIndex++;
     75//                                      }
     76//                                      count++;
     77//                              }
     78//                      }
    6279        }
    6380

  • issm/trunk-jpl/src/m/extrusion/project3d.js

    r20838 r22098  
    1 project3d = function() {
    2     //PROJECT3D - vertically project a vector from 2d mesh
    3     //
    4     //   vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.
    5     //   This vector can be a node vector of size (md.mesh.numberofvertices2d,N/A) or an
    6     //   element vector of size (md.mesh.numberofelements2d,N/A).
    7     //   arguments:
    8     //      'vector': 2d vector
    9     //      'type': 'element' or 'node'.
    10     //   options:
    11     //      'layer' a layer number where vector should keep its values. If not specified, all layers adopt the
    12     //             value of the 2d vector.
    13     //      'padding': default to 0 (value adopted by other 3d layers not being projected0
    14     //
    15     //   Egs:
    16     // md.extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',null);
    17     // md.extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0);
    18     // md.extruded_vector=project3d(md,'vector',vector2d,'type','node');
     1function project3d() {
     2        //PROJECT3D - vertically project a vector from 2d mesh
     3        //
     4        //   vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.
     5        //   This vector can be a node vector of size (md.mesh.numberofvertices2d,N/A) or an
     6        //   element vector of size (md.mesh.numberofelements2d,N/A).
     7        //   arguments:
     8        //      'vector': 2d vector
     9        //      'type': 'element' or 'node'.
     10        //   options:
     11        //      'layer' a layer number where vector should keep its values. If not specified, all layers adopt the
     12        //             value of the 2d vector.
     13        //      'padding': default to 0 (value adopted by other 3d layers not being projected0
     14        //
     15        //   Egs:
     16        // md.extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',null);
     17        // md.extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0);
     18        // md.extruded_vector=project3d(md,'vector',vector2d,'type','node');
    1919
    20     //some regular checks
    21    
    22     function remove_first_n(start, arglist) { // Using slice() on arguments is discouraged because it prevents optimizations in engines such as V8
    23         var args = [];
     20        //some regular checks
     21        if (arguments.length===1 || arguments.length===0) {
     22                console.error('project3d bad usage');
     23        }
     24        if (md.mesh.elementtype() !== 'Penta') {
     25                console.error('input model is not 3d');
     26        }
    2427
    25         for (var i = start; i < arglist.length; i++) {
    26             args.push(arglist[i]);
    27         }
     28        //retrieve parameters from options.
     29        var options      = new pairoptions(Array.prototype.slice.call(arguments, 1)); // slice to remove md
     30        var vector2d     = options.getfieldvalue('vector');     //mandatory
     31        var type         = options.getfieldvalue('type');       //mandatory
     32        var layer        = options.getfieldvalue('layer',0);    //optional (do all layers default:)
     33        var paddingvalue = options.getfieldvalue('padding',0);  //0 by default
    2834
    29         return args;
    30     }
     35        if (Number.isNaN(vector2d) || vector2d === 0 || vector2d.length === 1) { // NaN treated as length 1 in MATLAB
     36                projected_vector=vector2d;
     37        } else if (type.toLowerCase() === 'node') {
     38                //Initialize 3d vector
     39                if (vector2d.length===md.mesh.numberofvertices2d) {
     40                        projected_vector=NewArrayFill(md.mesh.numberofvertices,paddingvalue);
     41                } else if (vector2d.length===md.mesh.numberofvertices2d+1) {
     42                        projected_vector=NewArrayFill(md.mesh.numberofvertices+1,paddingvalue);
     43                        projected_vector[projected_vector.length-1] = vector2d[vector2d.length-1];
     44                        vector2d.pop();
     45                } else {
     46                        console.error('vector length not supported')
     47                }
    3148
    32     if (arguments.length===1 || arguments.length===0) {
    33         //help project3d
    34         console.error('project3d bad usage');
    35     }
     49                //Fill in
     50                if (layer===0) {
     51                        for (var i = 1; i <= md.mesh.numberoflayers; ++i) {
     52                                for (var j = (i-1)*md.mesh.numberofvertices2d, k = 0; j < i*md.mesh.numberofvertices2d; j++, k++) {
     53                                        projected_vector[j] = vector2d[k];
     54                                }
     55                        }
     56                } else {
     57                        for (var j = (layer-1)*md.mesh.numberofvertices2d, k = 0; j < layer*md.mesh.numberofvertices2d; j++, k++) {
     58                                projected_vector[j] = vector2d[k];
     59                        }
     60                }
     61        } else if (type.toLowerCase() === 'element') {
     62                //Initialize 3d vector
     63                //var vector2d_size2 = Array.isArray(vector2d[0]) ? vector2d[0].length : 1; //get size of vector2d's 2nd axis
     64                if (vector2d.length===md.mesh.numberofelements2d) {
     65                        projected_vector=NewArrayFill(md.mesh.numberofelements,paddingvalue);
     66                } else if (vector2d.length===md.mesh.numberofelements2d+1) {
     67                        projected_vector=NewArrayFill(md.mesh.numberofelements+1,paddingvalue);
     68                        projected_vector[projected_vector.length-1] = vector2d[vector2d.length-1];
     69                        vector2d.pop();
     70                } else {
     71                        console.error('vector length not supported')
     72                }
    3673
    37     var md = arguments[0];
     74                //Fill in
     75                if (layer===0) {
     76                        for (var i = 1; i <= md.mesh.numberoflayers-1; ++i) {
     77                                for (var j = (i-1)*md.mesh.numberofelements2d, k = 0; j < i*md.mesh.numberofelements2d; j++, k++) {
     78                                        projected_vector[j] = vector2d[k];
     79                                }
     80                        }
     81                } else {
     82                        for (var j = (layer-1)*md.mesh.numberofelements2d, k = 0; j < layer*md.mesh.numberofelements2d; j++, k++) {
     83                                projected_vector[j] = vector2d[k];
     84                        }
     85                }
     86        } else {
     87                console.error('project3d error message: unknown projection type');
     88        }
    3889
    39     if (md.mesh.elementtype() !== 'Penta') {
    40         console.error('input model is not 3d');
    41     }
    42 
    43     //retrieve parameters from options.
    44     options      = new pairoptions(remove_first_n(1, arguments)); // slice to remove md
    45     vector2d     = options.getfieldvalue('vector');     //mandatory
    46     type         = options.getfieldvalue('type');       //mandatory
    47     layer        = options.getfieldvalue('layer',0);    //optional (do all layers default:)
    48     paddingvalue = options.getfieldvalue('padding',0);  //0 by default
    49 
    50     if (isNaN(vector2d) || vector2d === 0 || vector2d.length === 1) { // NaN treated as length 1 in MATLAB
    51         projected_vector=vector2d;
    52     } else if (type.toLowerCase() === 'node') {
    53 
    54         //Initialize 3d vector
    55         if (vector2d.length===md.mesh.numberofvertices2d) {
    56                 projected_vector=ones(md.mesh.numberofvertices,  vector2d[0].length).map(function(arr) { return arr.fill(paddingvalue); });
    57         } else if (vector2d.length===md.mesh.numberofvertices2d+1) {
    58             projected_vector=ones(md.mesh.numberofvertices+1,vector2d[0].length).map(function(arr) { return arr.fill(paddingvalue); });
    59             projected_vector[projected_vector.length-1] = projected_vector[projected_vector.length-1].map(function(element, index) { return vector2d[vector2d.length-1][index]; });
    60             vector2d.pop(); // Remove last array
    61         } else {
    62             console.error('vector length not supported')
    63         }
    64 
    65             //Fill in
    66         if (layer===0) {
    67             for (var i = 1; i < md.mesh.numberoflayers; ++i) {
    68                 var vec_idx = 0;
    69 
    70                 for (var j = (i-1)*md.mesh.numberofvertices2d+1, vec_idx = 0; j <= i * md.mesh.numberofvertices2d && vec_idx < vector2d.length; ++j, ++vec_idx) {
    71                     projected_vector[j].map(function(element, index) { return vector2d[vec_idx][index]; });
    72                 }
    73 
    74             }
    75         } else {
    76             var vec_idx = 0;
    77 
    78             for (var i = (layer-1)*md.mesh.numberofvertices2d+1, vec_idx = 0; i <= layer * md.mesh.numberofvertices2d && vec_idx < vector2d.length; ++i, ++vec_idx) {
    79                 projected_vector[i] = vector2d[vec_idx];
    80             }
    81         }
    82     } else if (type.toLowerCase() === 'element') {
    83             //Initialize 3d vector
    84         if (vector2d.length===md.mesh.numberofelements2d) {
    85                 projected_vector = ones(md.mesh.numberofelements,  vector2d[0].length).map(function(arr) { return arr.fill(paddingvalue); });
    86         } else if (vector2d.length===md.mesh.numberofelements2d+1) {
    87             projected_vector = ones(md.mesh.numberofelements+1,  vector2d[0].length).map(function(arr) { return arr.fill(paddingvalue); });
    88 
    89             projected_vector[projected_vector.length-1].map(function(element, index) { return vector2d[vector2d.length-1][index]; });
    90 
    91             vector2d.pop();
    92         } else {
    93             console.error('vector length not supported')
    94         }
    95 
    96             //Fill in
    97         if (layer===0) {
    98             for (var i = 1; i < md.mesh.numberoflayers-1; ++i) {
    99                 var vec_idx=0;
    100                 for (var j = (i-1)*md.mesh.numberofelements2d+1, vec_idx = 0; j <= i * md.mesh.numberofelements2d && vec_idx < vector2d.length; ++j, ++vec_idx) {
    101                     projected_vector[j].map(function(element, index) { return vector2d[vec_idx][index]; });
    102                 }
    103             }
    104         } else {
    105             var vec_idx=0;
    106             for (var i = (layer-1)*md.mesh.numberofelements2d+1, vec_idx = 0; i <= layer * md.mesh.numberofelements2d && vec_idx < vector2d.length; ++i, ++vec_idx) {
    107                 projected_vector[i].map(function(element, index) { return vector2d[vec_idx][index]; });
    108             }
    109         }
    110     } else {
    111         console.error('project3d error message: unknown projection type');
    112     }
    113 
    114     return projected_vector;
     90        return projected_vector;
    11591};
Note: See TracChangeset for help on using the changeset viewer.