Changeset 26308


Ignore:
Timestamp:
06/08/21 16:38:39 (4 years ago)
Author:
jdquinn
Message:

CHG: MATLAB -> JS

Location:
issm/trunk-jpl/src/m
Files:
8 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r24469 r26308  
    1 //dsl Class definition
    2 //
    3 //   Usage:
    4 //      dsl=dsl();
     1class dsl {//{{{
     2        /**
     3         * DSL class definition
     4         *
     5         * Usage:
     6         *     dsl=dsl();
     7         *
     8         * TODO:
     9         * - Implement from struct constructor (see dsl.m)?
     10         */
     11        constructor() {//{{{
     12                this.global_average_thermosteric_sea_level; // Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m)
     13                this.sea_surface_height_above_geoid; // Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m)
     14                this.sea_water_pressure_at_sea_floor; // Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!)
    515
    6 function dsl(){
    7         //methods
    8         this.setdefaultparameters = function(){// {{{
    9                 this.requested_outputs=['default'];
    10         } // }}}
    11         this.disp = function(){ // {{{
    12                 console.log(sprintf('   dsl parameters:'));
    13                 fielddisplay(this,'global_average_thermosteric_sea_level_change','corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable global rate (mm/yr)');
    14                 fielddisplay(this,'sea_surface_height_change_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)');
    15                 fielddisplay(this,'sea_water_pressure_change_at_sea_floor','corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)');
     16                let nargs = arguments.length;
     17                if (nargs == 0) {
     18                        this.setdefaultparameters();
     19                } else {
     20                        error('constructor not supported');
     21                }
     22        } //}}}
    1623
     24        setdefaultparameters() {//{{{
     25                this.global_average_thermosteric_sea_level = NaN;
     26                this.sea_surface_height_above_geoid = NaN;
     27                this.sea_water_pressure_at_sea_floor = NaN;
     28        } //}}}
    1729
    18         } // }}}
    19         this.defaultoutputs = function(){ // {{{
    20                 return '';
    21         }//}}}
    22     this.classname = function(){ // {{{
    23         return "dsl";
    24     } // }}}
    25     this.extrude = function(md) {//{{{
    26         this.sea_surface_height_change_above_geoid=project3d(md,'vector',this.sea_surface_height_change_above_geoid,'type','node');
    27         this.sea_water_pressure_change_at_sea_floor=project3d(md,'vector',this.sea_water_pressure_change_at_sea_floor,'type','node');
    28         return this;
    29     }//}}}
    30     this.checkconsistency = function(md,solution,analyses) { //{{{
     30        checkconsistency(md, solution, analyses) {//{{{
     31                if (!analyses.includes('SealevelchangeAnalysis') || (solution === 'TransientSolution' && !md.transient.isslc) || !md.transient.isoceantransport) {
     32                        return md;
     33                }
    3134
    32         if(ArrayAnyEqual(ArrayIsMember('SealevelriseAnalysis',analyses),1)){
    33             checkfield(md,'fieldname','dsl.sea_surface_height_change_above_geoid','timeseries',1,'NaN',1,'Inf',1);
    34             checkfield(md,'fieldname','dsl.sea_water_pressure_change_at_sea_floor','timeseries',1,'NaN',1,'Inf',1);
    35         }
     35                md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1);
     36                md = checkfield(md, 'fieldname', 'dsl.sea_surface_height_above_geoid', 'NaN', 1, 'Inf', 1, 'timeseries', 1);
     37                md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'NaN', 1, 'Inf', 1, 'timeseries', 1);
    3638
    37     } // }}}
    38     this.marshall=function(md,prefix,fid) { //{{{
     39                if (md.solidearth.settings.compute_bp_grd) {
     40                        md = checkfield(md, 'fieldname', 'dsl.sea_water_pressure_at_sea_floor', 'empty', 1);
     41                }
    3942
    40         var yts=md.constants.yts;
     43                return md;
     44        } //}}}
    4145
    42                 WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer');
    43                 WriteData(fid,prefix,'object',this,'class','dsl','fieldname','global_average_thermosteric_sea_level_change','format','DoubleMat','mattype',1,'timeserieslength',1+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts);
    44                         WriteData(fid,prefix,'object',this,'class','dsl','fieldname','sea_water_pressure_change_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts);
    45                         WriteData(fid,prefix,'object',this,'class',dsl,'fieldname','sea_surface_height_change_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    46     }//}}}
    47     this.fix=function() { //{{{
    48     }//}}}
    49         //properties
    50     // {{{
    51         this.global_average_thermosteric_sea_level_change = NaN;
    52         this.sea_surface_height_change_above_geoid = NaN;
    53         this.sea_water_pressure_change_at_sea_floor = NaN;
    54     // }}}
    55 }
     46        disp() {//{{{
     47                console.log('WARNING: dsl::disp is not yet implemented');
     48        } //}}}
     49
     50        marshall(md, prefix, fid) {//{{{
     51                WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer');
     52                WriteData(fid, prefix, 'object', this, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeseries', 1, 'timeserieslength', 2, 'yts', md.constants.yts); // mattype 2, because we are sending a GMSL value identical everywhere on each element
     53                WriteData(fid, prefix, 'object', this, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts); // mattype 1 because we specify DSL at vertex locations
     54                WriteData(fid, prefix, 'object', this, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts); // mattype 1 because we specify bottom pressure at vertex locations
     55        } //}}}
     56
     57        extrude(md) {//{{{
     58                this.sea_surface_height_above_geoid = project3d(md, 'vector', this.sea_surface_height_above_geoid, 'type', 'node', 'layer', 1);
     59                this.sea_water_pressure_at_sea_floor = project3d(md, 'vector', this.sea_water_pressure_at_sea_floor, 'type', 'node', 'layer', 1);
     60        } //}}}
     61} //}}}

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

    r22283 r26308  
    1 //HYDROLOGYSHREVE class definition
    2 //
    3 //   Usage:
    4 //      hydrologyshreve=new hydrologyshreve();
     1class hydrologyshreve {//{{{
     2        /**
     3         * HYDROLOGYSHREVE class definition
     4         *
     5         * Usage:
     6         *     hydrologyshreve = hydrologyshreve();
     7         */
     8        constructor() {//{{{
     9                this.spcwatercolumn = NaN;
     10                this.stabilization = 0;
     11                this.requested_outputs = [];
    512
    6 function hydrologyshreve (){
    7         //methods
    8         this.setdefaultparameters = function(){// {{{
     13                if (arguments.length == 0) {
     14                        this.setdefaultparameters();
     15                } else {
     16                        error('constructor not supported');
     17                }
     18        } //}}}
    919
    10                 //Type of stabilization to use 0:nothing 1:artificial_diffusivity
    11                 this.stabilization=1;
     20        disp() {//{{{
     21                console.log(sprintf('   hydrologyshreve solution parameters:'));
     22                fielddisplay(this, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]');
     23                fielddisplay(this, 'stabilization', 'artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.');
     24                fielddisplay(this, 'requested_outputs', 'additional outputs requested');
     25        } //}}}
    1226
    13                 //default output
    14                 this.requested_outputs=['default'];
     27        setdefaultparameters() {//{{{
     28                // Type of stabilization to use 0:nothing 1:artificial_diffusivity
     29                this.stabilization = 1;
     30                this.requested_outputs = ['default'];
     31        } //}}}
    1532
    16         }// }}}
    17                 this.disp= function(){// {{{
     33        checkconsistency(md, solution, analyses) {//{{{
     34                // Early return
     35                if (!analyses.contains('HydrologyShreveAnalysis') || (solution == 'TransientSolution' && !md.transient.ishydrology)) {
     36                        return md;
     37                }
     38                md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1);
     39                md = checkfield(md, 'fieldname', 'hydrology.stabilization', '>=', 0);
     40                return md;
     41        } //}}}
    1842
    19                 console.log(sprintf('   hydrologyshreve solution parameters:'));
    20                 fielddisplay(this,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
    21                 fielddisplay(this,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.');
    22                 fielddisplay(this,'requested_outputs','additional outputs requested');
     43        defaultoutputs(md) {//{{{
     44                return ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy'];
     45        } //}}}
    2346
    24         }// }}}
    25     this.extrude = function(md) {//{{{
    26         return this;
    27     };//}}}
    28                 this.classname= function(){// {{{
    29                         return "hydrologyshreve";
     47        marshall(md, prefix, fid) {//{{{
     48                WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer');
     49                WriteData(fid, prefix, 'object', this, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts);
     50                WriteData(fid, prefix, 'object', this, 'fieldname', 'stabilization', 'format', 'Double');
    3051
    31                 }// }}}
    32         this.type= function(){// {{{
     52                // Process requested outputs
     53                let outputs = this.requested_outputs;
     54                let pos = find(ismember(outputs, 'default'));
     55                if (pos.length) {
     56                        /*
     57                        NOTE: In order to handle case where user has added 'default' more
     58                        than once to this.requested_outputs, need to remove elements by
     59                        index in reverse order.
     60                        */
     61                        for (let i = (pos.length - 1); i >= 0; --i) {
     62                                outputs.splice(pos[i], 1);
     63                        }
     64                        let default_outputs = this.defaultoutputs(md);
     65                        for (let i = 0; i < default_outputs.length; ++i) {
     66                                outputs.push(default_outputs[i]);
     67                        }
     68                }
     69                WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray');
     70        } //}}}
    3371
    34                 return "hydrologyshreve";
    35         }// }}}
    36                 this.checkconsistency = function(md,solution,analyses) { //{{{
    37 
    38                         //Early return
    39                         if(!ArrayAnyEqual(ArrayIsMember('HydrologyShreveAnalysis',analyses),1)) return;
    40 
    41                         checkfield(md,'fieldname','hydrology.spcwatercolumn','Inf',1,'timeseries',1);
    42                         checkfield(md,'fieldname','hydrology.stabilization','>=',0);
    43                         checkfield(md,'fieldname','hydrology.requested_outputs','stringrow',1);
    44 
    45                 } // }}}
    46                 this.marshall=function(md,prefix,fid) { //{{{
    47                         WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
    48                         WriteData(fid,prefix,'object',this,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    49                         WriteData(fid,prefix,'object',this,'fieldname','stabilization','format','Double');
    50                         var outputs = this.requested_outputs;
    51                         for (var i=0;i<outputs.length;i++){
    52                                 if (outputs[i] == 'default') {
    53                                         outputs.splice(i,1);
    54                                         var newoutputs=this.defaultoutputs(md);
    55                                         for (var j=0;j<newoutputs.length;j++) outputs.push(newoutputs[j]);
    56                                 }
    57                         }
    58                         WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
    59                 }//}}}
    60                 this.defaultoutputs = function(md) { //{{{
    61                                 return ['Watercolumn','HydrologyWaterVx','HydrologyWaterVy'];
    62                 }//}}}
    63                 this.fix=function() { //{{{
    64                         this.spcwatercolumn=NullFix(this.spcwatercolumn,NaN);
    65                 }//}}}
    66         //properties
    67         // {{{
    68         this.spcwatercolumn = NaN;
    69         this.stabilization  = 0;
    70         this.requested_outputs = [];
    71         this.setdefaultparameters();
    72         //}}}
    73 }
     72        extrude(md) {//{{{
     73        } //}}}
     74} //}}}

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

    r25870 r26308  
    1 //MASSTRANSPORT class definition
    2 //
    3 //   Usage:
    4 //      masstransport=new masstransport();
     1class masstransport {//{{{
     2        /**
     3         * MASSTRANSPORT class definition
     4         *
     5         * Usage:
     6         *     masstransport = masstransport();
     7         */
     8        constructor() {//{{{
     9                this.spcthickness = NaN;
     10                this.isfreesurface = 0;
     11                this.min_thickness = 0;
     12                this.hydrostatic_adjustment = 0;
     13                this.stabilization = 0;
     14                this.vertex_pairing = 0;
     15                this.penalty_factor = 0;
     16                this.requested_outputs = 0;
    517
    6 function masstransport (){
    7         //methods
    8         this.setdefaultparameters = function(){// {{{
     18                if (arguments.length == 0) {
     19                        this.setdefaultparameters();
     20                } else {
     21                        error('constructor not supported');
     22                }
     23        } //}}}
    924
    10                 //Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin
    11                 this.stabilization=1;
     25        disp() {//{{{
     26                console.log('WARNING: masstransport::disp is not yet implemented');
     27        } //}}}
    1228
    13                 //Factor applied to compute the penalties kappa=max(stiffness matrix)*10^penalty_factor
    14                 this.penalty_factor=3;
     29        setdefaultparameters() {//{{{
     30                // Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin
     31                this.stabilization = 1;
    1532
    16                 //Minimum ice thickness that can be used
    17                 this.min_thickness=1;
     33                // Factor applied to compute the penalties kappa=max(stiffness matrix)*10^penalty_factor
     34                this.penalty_factor = 3;
    1835
    19                 //Hydrostatic adjustment
    20                 this.hydrostatic_adjustment='Absolute';
     36                // Minimum ice thickness that can be used
     37                this.min_thickness = 1;
    2138
    22                 //default output
    23                 this.requested_outputs=['default'];
     39                // Hydrostatic adjustment
     40                this.hydrostatic_adjustment = 'Absolute';
    2441
    25         }// }}}
    26         this.disp= function(){// {{{
    27                 console.log(sprintf('   Masstransport solution parameters:'));
    28                 fielddisplay(this,'spcthickness','thickness constraints (NaN means no constraint) [m]');
    29                 fielddisplay(this,'isfreesurface','do we use free surfaces (FS only) or mass conservation');
    30                 fielddisplay(this,'min_thickness','minimum ice thickness allowed [m]');
    31                 fielddisplay(this,'hydrostatic_adjustment',"adjustment of ice shelves surface and bed elevations: 'Incremental' or 'Absolute' ");
    32                 fielddisplay(this,'stabilization','0: no stabilization, 1: artificial diffusion, 2: streamline upwinding, 3: discontinuous Galerkin, 4: flux corrected transport, 5: streamline upwind Petrov-Galerkin (SUPG)');
     42                // Default output
     43                this.requested_outputs = ['default'];
     44        } //}}}
    3345
    34                 console.log(sprintf('\n      %s','Penalty options:'));
    35                 fielddisplay(this,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
    36                 fielddisplay(this,'vertex_pairing','pairs of vertices that are penalized');
    37                 fielddisplay(this,'requested_outputs','additional outputs requested');
     46        checkconsistency(md, solution, analyses) {//{{{
     47                // Early return
     48                if (analyses.includes('HydrologyShreveAnalysis') || (solution == 'TransientSolution' && !md.trans.ismasstransport)) {
     49                        return md;
     50                }
    3851
    39         }// }}}
    40         this.classname= function(){// {{{
    41                 return "masstransport";
    42         }// }}}
    43     this.extrude = function(md) {//{{{
    44         this.spcthickness=project3d(md,'vector',this.spcthickness,'type','node');
    45         return this;
    46     }//}}}
    47                 this.checkconsistency = function (md,solution,analyses){  // {{{
     52                md = checkfield(md, 'fieldname', 'masstransport.spcthickness', 'Inf', 1, 'timeseries', 1);
     53                md = checkfield(md, 'fieldname', 'masstransport.isfreesurface', 'values', [0, 1]);
     54                md = checkfield(md, 'fieldname', 'masstransport.hydrostatic_adjustment', 'values', ['Absolute', 'Incremental']);
     55                md = checkfield(md, 'fieldname', 'masstransport.stabilization', 'values', [0, 1, 2, 3, 4, 5]);
     56                md = checkfield(md, 'fieldname', 'masstransport.min_thickness', '>', 0);
     57                md = checkfield(md, 'fieldname', 'masstransport.requested_outputs', 'stringrow', 1);
     58                if (!any(isnan(md.stressbalance.vertex.vertex_pairing))) {
     59                        md = checkfield(md, 'fieldname', 'stressbalance.vertex_pairing', '>', 0);
     60                }
    4861
    49                         //Early return
    50                         if(!ArrayAnyEqual(ArrayIsMember('HydrologyShreveAnalysis',analyses),1) | (solution=='TransientSolution' & md.trans.ismasstransport==0)) return;
     62                return md;
     63        } //}}}
    5164
    52                         checkfield(md,'fieldname','masstransport.spcthickness','Inf',1,'timeseries',1);
    53                         checkfield(md,'fieldname','masstransport.isfreesurface','values',[0 ,1]);
    54                         checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',['Absolute', 'Incremental']);
    55                         checkfield(md,'fieldname','masstransport.stabilization','values',[0,1,2,3,4]);
    56                         checkfield(md,'fieldname','masstransport.min_thickness','>',0);
    57                         checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1);
     65        defaultoutputs(md) {//{{{
     66                return ['Thickness', 'Surface', 'Base'];
     67        } //}}}
    5868
    59                 } // }}}
    60                 this.marshall=function(md,prefix,fid) { //{{{
     69        marshall(md, prefix, fid) {//{{{
     70                let yts = md.constants.yts;
    6171
    62                         var yts=md.constants.yts;
     72                WriteData(fid, prefix, 'object', this, 'fieldname', 'spcthickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
     73                WriteData(fid, prefix, 'object', this, 'fieldname', 'isfreesurface', 'format', 'Boolean');
     74                WriteData(fid, prefix, 'object', this, 'fieldname', 'min_thickness', 'format', 'Double');
     75                WriteData(fid, prefix, 'data', this.hydrostatic_adjustment, 'format', 'String', 'name', 'md.masstransport.hydrostatic_adjustment');
     76                WriteData(fid, prefix, 'object', this, 'fieldname', 'stabilization', 'format', 'Integer');
     77                WriteData(fid, prefix, 'object', this, 'fieldname', 'vertex_pairing', 'format', 'DoubleMat', 'mattype', 3);
     78                WriteData(fid, prefix, 'object', this, 'fieldname', 'penalty_factor', 'format', 'Double');
    6379
    64                         WriteData(fid,prefix,'object',this,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    65                         WriteData(fid,prefix,'object',this,'fieldname','isfreesurface','format','Boolean');
    66                         WriteData(fid,prefix,'object',this,'fieldname','min_thickness','format','Double');
    67                         WriteData(fid,prefix,'data',this.hydrostatic_adjustment,'format','String','name','md.masstransport.hydrostatic_adjustment');
    68                         WriteData(fid,prefix,'object',this,'fieldname','stabilization','format','Integer');
    69                         WriteData(fid,prefix,'object',this,'fieldname','vertex_pairing','format','DoubleMat','mattype',3);
    70                         WriteData(fid,prefix,'object',this,'fieldname','penalty_factor','format','Double');
     80                // Process requested outputs
     81                let outputs = this.requested_outputs;
     82                let pos = find(ismember(outputs, 'default'));
     83                if (pos.length) {
     84                        /*
     85                        NOTE: In order to handle case where user has added 'default' more
     86                        than once to this.requested_outputs, need to remove elements by
     87                        index in reverse order.
     88                        */
     89                        for (let i = (pos.length - 1); i >= 0; --i) {
     90                                outputs.splice(pos[i], 1);
     91                        }
     92                        let default_outputs = this.defaultoutputs(md);
     93                        for (let i = 0; i < default_outputs.length; ++i) {
     94                                outputs.push(default_outputs[i]);
     95                        }
     96                }
     97                WriteData(fid, prefix, 'data', outputs, 'name', 'md.masstransport.requested_outputs', 'format', 'StringArray');
     98        } //}}}
    7199
    72                         //process requested outputs
    73                         var outputs = this.requested_outputs;
    74                         for (var i=0;i<outputs.length;i++){
    75                                 if (outputs[i] == 'default') {
    76                                         outputs.splice(i,1);
    77                                         var newoutputs=this.defaultoutputs(md);
    78                                         for (var j=0;j<newoutputs.length;j++) outputs.push(newoutputs[j]);
    79                                 }
    80                         }
    81                         WriteData(fid,prefix,'data',outputs,'name','md.masstransport.requested_outputs','format','StringArray');
    82                 }//}}}
    83                 this.defaultoutputs = function(md) { //{{{
    84                         return ['Thickness','Surface','Base'];
    85                 }//}}}
    86                 this.fix=function() { //{{{
    87                         this.spcthickness=NullFix(this.spcthickness,NaN);
    88                         this.vertex_pairing=NullFix(this.vertex_pairing,NaN);
    89                 }//}}}
    90         //properties
    91         // {{{
    92 
    93         this.spcthickness           = NaN;
    94         this.isfreesurface          = 0;
    95         this.min_thickness          = 0;
    96         this.hydrostatic_adjustment = 0;
    97         this.stabilization          = 0;
    98         this.vertex_pairing         = NaN;
    99         this.penalty_factor         = 0;
    100         this.requested_outputs      = [];
    101 
    102         this.setdefaultparameters();
    103         //}}}
     100        extrude(md) {//{{{
     101        } //}}}
    104102}

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

    r26306 r26308  
    1 //MATERIALS class definition
    2 //
    3 //      Usage:
    4 //              materials=materials();
    5 
    6 class materials {
    7         //methods
     1class materials {//{{{
     2        /**
     3         * MATERIALS class definition
     4         *
     5         * Usage:
     6         *     materials = materials();
     7         */
    88        constructor() {//{{{
    99                this.nature = [];
    … …  
    1111                let nargs = arguments.length;
    1212                if (nargs == 0) {
    13                         this.nature=['ice'];
     13                        this.nature = ['ice'];
    1414                } else {
    15                         this.nature=arguments;
    16                 }
    17 
    18                 //check this is acceptable:
    19                 for (let i=0; i < length(this.nature); i++){
    20                         if (!(strcmpi(this.nature[i],'litho') || strcmpi(this.nature[i],'ice') || strcmpi(this.nature[i],'hydro'))) {
    21                                 error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
    22                         }
    23                 }
    24 
    25                 //start filling in the dynamic fields:
    26                 for (let i=0; i < length(this.nature); i++) {
    27                         let nat=this.nature[i];
    28                         switch (nat) {
    29                                 case 'ice':
    30                                         this.rho_ice = 0;
    31                                         this.rho_water = 0;
    32                                         this.rho_freshwater = 0;
    33                                         this.mu_water = 0;
    34                                         this.heatcapacity = 0;
    35                                         this.latentheat = 0;
    36                                         this.thermalconductivity = 0;
    37                                         this.temperateiceconductivity = 0;
    38                                         this.effectiveconductivity_averaging = 0;
    39                                         this.meltingpoint = 0;
    40                                         this.beta = 0;
    41                                         this.mixed_layer_capacity = 0;
    42                                         this.thermal_exchange_velocity = 0;
    43                                         this.rheology_B = 0;
    44                                         this.rheology_n = 0;
    45                                         this.rheology_law = 0;
    46                                         break;
    47                                 case 'litho':
    48                                         this.numlayers = 0;
    49                                         this.radius = 0;
    50                                         this.viscosity = 0;
    51                                         this.lame_lambda = 0;
    52                                         this.lame_mu = 0;
    53                                         this.burgers_viscosity = 0;
    54                                         this.burgers_mu = 0;
    55                                         this.ebm_alpha = 0;
    56                                         this.ebm_delta = 0;
    57                                         this.ebm_taul = 0;
    58                                         this.ebm_tauh = 0;
    59                                         this.rheologymodel = 0;
    60                                         this.density = 0;
    61                                         this.issolid = 0;
    62                                         break;
    63                                 case 'hydro':
    64                                         this.rho_ice = 0;
    65                                         this.rho_water = 0;
    66                                         this.rho_freshwater = 0;
    67                                         break;
    68                                 default:       
    69                                         error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
     15                        this.nature = arguments;
     16                }
     17
     18                // Check this is acceptable
     19                for (let i = 0; i < this.nature.length; ++i) {
     20                        if (!(strcmpi(this.nature[i], 'litho') || strcmpi(this.nature[i], 'ice') || strcmpi(this.nature[i], 'hydro'))) {
     21                                error('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')');
     22                        }
     23                }
     24
     25                // Start filling in the dynamic fields (not truly dynamic under JavaScript)
     26                for (let i = 0; i < this.nature.length; ++i) {
     27                        let nat = this.nature[i];
     28                        if (nat == 'ice') {
     29                                this.rho_ice = 0;
     30                                this.rho_water = 0;
     31                                this.rho_freshwater = 0;
     32                                this.mu_water = 0;
     33                                this.heatcapacity = 0;
     34                                this.latentheat = 0;
     35                                this.thermalconductivity = 0;
     36                                this.temperateiceconductivity = 0;
     37                                this.effectiveconductivity_averaging = 0;
     38                                this.meltingpoint = 0;
     39                                this.beta = 0;
     40                                this.mixed_layer_capacity = 0;
     41                                this.thermal_exchange_velocity = 0;
     42                                this.rheology_B = 0;
     43                                this.rheology_n = 0;
     44                                this.rheology_law = 0;
     45                        } else if (nat == 'litho') {
     46                                this.numlayers = 0;
     47                                this.radius = 0;
     48                                this.viscosity = 0;
     49                                this.lame_lambda = 0;
     50                                this.lame_mu = 0;
     51                                this.burgers_viscosity = 0;
     52                                this.burgers_mu = 0;
     53                                this.ebm_alpha = 0;
     54                                this.ebm_delta = 0;
     55                                this.ebm_taul = 0;
     56                                this.ebm_tauh = 0;
     57                                this.rheologymodel = 0;
     58                                this.density = 0;
     59                                this.issolid = 0;
     60                        } else if (nat == 'hydro') {
     61                                this.rho_ice = 0;
     62                                this.rho_water = 0;
     63                                this.rho_freshwater = 0;
     64                        } else {
     65                                error('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')');
    7066                        }
    7167                }
    7268                this.earth_density = 0;
    7369
    74                 //set default parameters:
     70                // Set default parameters
    7571                this.setdefaultparameters();
    7672        } //}}}
     73
     74        disp() {//{{{
     75                console.log('WARNING: materials::disp is not yet implemented');
     76        } //}}}
     77
    7778        setdefaultparameters() {//{{{
    78                 for (let i=0; i < length(this.nature); i++) {
    79                         let nat=this.nature[i];
    80                         switch (nat) {
    81                                 case 'ice':
    82                                         //ice density (kg/m^3)
    83                                         this.rho_ice=917.;
    84 
    85                                         //ocean water density (kg/m^3)
    86                                         this.rho_water=1023.;
    87 
    88                                         //fresh water density (kg/m^3)
    89                                         this.rho_freshwater=1000.;
    90 
    91                                         //water viscosity (N.s/m^2)
    92                                         this.mu_water=0.001787;
    93 
    94                                         //ice heat capacity cp (J/kg/K)
    95                                         this.heatcapacity=2093.;
    96 
    97                                         //ice latent heat of fusion L (J/kg)
    98                                         this.latentheat=3.34*1e5;
    99 
    100                                         //ice thermal conductivity (W/m/K)
    101                                         this.thermalconductivity=2.4;
    102 
    103                                         //wet ice thermal conductivity (W/m/K)
    104                                         this.temperateiceconductivity=.24;
    105 
    106                                         //computation of effective conductivity
    107                                         this.effectiveconductivity_averaging=1;
    108 
    109                                         //the melting point of ice at 1 atmosphere of pressure in K
    110                                         this.meltingpoint=273.15;
    111 
    112                                         //rate of change of melting point with pressure (K/Pa)
    113                                         this.beta=9.8*1e-8;
    114 
    115                                         //mixed layer (ice-water interface) heat capacity (J/kg/K)
    116                                         this.mixed_layer_capacity=3974.;
    117 
    118                                         //thermal exchange velocity (ice-water interface) (m/s)
    119                                         this.thermal_exchange_velocity=1.00*1e-4;
    120 
    121                                         //Rheology law: what is the temperature dependence of B with T
    122                                         //available: none, paterson and arrhenius
    123                                         this.rheology_law='Paterson';
    124 
    125                                         //Rheology fields default:
    126                                         this.rheology_B   = 1 * 1e8;
    127                                         this.rheology_n   = 3;
    128                     break;
    129 
    130                                 case 'litho':
    131                                         //we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
    132                                         this.numlayers=2;
    133 
    134                                         //center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface
    135                                         //(with 1d3 to avoid numerical singularities)
    136                                         this.radius=[1e3,6278*1e3,6378*1e3];
    137 
    138                                         this.viscosity=[1e21,1e40]; //mantle and lithosphere viscosity (respectively) [Pa.s]
    139                                         this.lame_mu=[1.45*1e11,6.7*1e10];  // (Pa) //lithosphere and mantle shear modulus (respectively) [Pa]
    140                                         this.lame_lambda=this.lame_mu;  // (Pa) //mantle and lithosphere lamba parameter (respectively) [Pa]
    141                                         this.burgers_viscosity=[NaN,NaN];
    142                                         this.burgers_mu=[NaN,NaN];
    143 
    144                                         this.ebm_alpha=[NaN,NaN];
    145                                         this.ebm_delta=[NaN,NaN];
    146                                         this.ebm_taul=[NaN,NaN];
    147                                         this.ebm_tauh=[NaN,NaN];
    148                                         this.rheologymodel=[0,0];
    149                                         this.density=[5.51*1e3,5.50*1e3];  // (Pa) //mantle and lithosphere density [kg/m^3]
    150                                         this.issolid=[1,1]; // is layer solid or liquid.
    151                     break;
    152 
    153                                 case 'hydro':
    154                                         //ice density (kg/m^3)
    155                                         this.rho_ice=917.;
    156 
    157                                         //ocean water density (kg/m^3)
    158                                         this.rho_water=1023.;
    159                                        
    160                                         //fresh water density (kg/m^3)
    161                                         this.rho_freshwater=1000.;
    162                     break;
    163 
    164                                 default:
    165                                         error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
    166                         }
    167 
    168                         // average density of the Earth (kg/m^3)
     79                for (let i = 0; i < this.nature.length; ++i) {
     80                        let nat = this.nature[i];
     81                        if (nat == 'ice') {
     82                                // Ice density (kg/m^3)
     83                                this.rho_ice = 917;
     84
     85                                // Ocean water density (kg/m^3)
     86                                this.rho_water = 1023
     87
     88                                // Fresh water density (kg/m^3)
     89                                this.rho_freshwater = 1000;
     90
     91                                // Water viscosity (N.s/m^2)
     92                                this.mu_water = 0.001787;
     93
     94                                // Ice heat capacity cp (J/kg/K)
     95                                this.heatcapacity = 2093;
     96
     97                                // Ice latent heat of fusion L (J/kg)
     98                                this.latentheat = 3.34 * 1e5;
     99
     100                                // Ice thermal conductivity (W/m/K)
     101                                this.thermalconductivity = 2.4;
     102
     103                                // Wet ice thermal conductivity (W/m/K)
     104                                this.temperateiceconductivity = 0.24;
     105
     106                                // Computation of effective conductivity
     107                                this.effectiveconductivity_averaging = 1;
     108
     109                                // The melting point of ice at 1 atmosphere of pressure in K
     110                                this.meltingpoint = 273.15;
     111
     112                                // Rate of change of melting point with pressure (K/Pa)
     113                                this.beta = 9.8 * 1e-8;
     114
     115                                // Mixed layer (ice-water interface) heat capacity (J/kg/K)
     116                                this.mixed_layer_capacity = 3974;
     117
     118                                // Thermal exchange velocity (ice-water interface) (m/s)
     119                                this.thermal_exchange_velocity = 1.00 * 1e-4;
     120
     121                                // Rheology law: what is the temperature dependence of B with T
     122                                // available: none, paterson and arrhenius
     123                                this.rheology_law = 'Paterson';
     124
     125                                // Rheology fields default
     126                                this.rheology_B = 1 * 1e8;
     127                                this.rheology_n = 3;
     128                        } else if (nat == 'litho') {
     129                                // We default to a configuration that enables running GIA
     130                                // solutions using giacaron and/or giaivins
     131                                this.numlayers = 2;
     132
     133                                // Center of the earth (approximation, must not be 0), then the
     134                                // lab (lithosphere/asthenosphere boundary) then the surface
     135                                // (with 1d3 to avoid numerical singularities)
     136                                this.radius = [1e3, 6278 * 1e3, 6378 * 1e3];
     137
     138                                this.viscosity = [1e21, 1e40]; // Mantle and lithosphere viscosity (respectively) [Pa.s]
     139                                this.lame_mu = [1.45 * 1e11, 6.7 * 1e10]; // (Pa) // Lithosphere and mantle shear modulus (respectively) [Pa]
     140                                this.lame_lambda = this.lame_mu; // (Pa) // Mantle and lithosphere lamba parameter (respectively) [Pa]
     141                                this.burgers_viscosity = [NaN, NaN];
     142                                this.burgers_mu = [NaN, NaN];
     143
     144                                this.ebm_alpha = [NaN, NaN];
     145                                this.ebm_delta = [NaN, NaN];
     146                                this.ebm_taul = [NaN, NaN];
     147                                this.ebm_tauh = [NaN, NaN];
     148                                this.rheologymodel = [0, 0];
     149                                this.density = [5.51 * 1e3, 5.50 * 1e3]; // (Pa) // Mantle and lithosphere density [kg/m^3]
     150                                this.issolid = [1, 1]; // Is layer solid or liquid?
     151                        } else if (nat == 'hydro') {
     152                                // Ice density (kg/m^3)
     153                                this.rho_ice = 917;
     154
     155                                // Ocean water density (kg/m^3)
     156                                this.rho_water = 1023;
     157
     158                                // Fresh water density (kg/m^3)
     159                                this.rho_freshwater = 1000;
     160                        } else {
     161                                error('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')');
     162                        }
     163
     164                        // Average density of the Earth (kg/m^3)
    169165                        this.earth_density=5512;
    170 
    171                 }
    172         } //}}}
    173         disp() {//{{{
    174                 console.log(sprintf('   Materials:'));
    175 
    176                 for (let i=0; i < length(this.nature); i++) {
    177                         let nat=this.nature[i];
    178                         switch (nat) {
    179                                 case 'ice':
    180                                         console.log(sprintf('   \nIce:'));
    181                                         fielddisplay(this,'rho_ice','ice density [kg/m^3]');
    182                                         fielddisplay(this,'rho_water','ocean water density [kg/m^3]');
    183                                         fielddisplay(this,'rho_freshwater','fresh water density [kg/m^3]');
    184                                         fielddisplay(this,'mu_water','water viscosity [N s/m^2]');
    185                                         fielddisplay(this,'heatcapacity','heat capacity [J/kg/K]');
    186                                         fielddisplay(this,'thermalconductivity','ice thermal conductivity [W/m/K]');
    187                                         fielddisplay(this,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
    188                                         fielddisplay(this,'meltingpoint','melting point of ice at 1atm in K');
    189                                         fielddisplay(this,'latentheat','latent heat of fusion [J/kg]');
    190                                         fielddisplay(this,'beta','rate of change of melting point with pressure [K/Pa]');
    191                                         fielddisplay(this,'mixed_layer_capacity','mixed layer capacity [W/kg/K]');
    192                                         fielddisplay(this,'thermal_exchange_velocity','thermal exchange velocity [m/s]');
    193                                         fielddisplay(this,'rheology_B','flow law parameter [Pa s^(1/n)]');
    194                                         fielddisplay(this,'rheology_n','Glen\'s flow law exponent');
    195                                         fielddisplay(this,'rheology_law','law for the temperature dependance of the rheology: "None", "BuddJacka", Cuffey", "CuffeyTemperate", "Paterson", "Arrhenius", "LliboutryDuval", "NyeCO2", or "NyeH2O"');
    196                                     break;
    197                                 case 'litho':
    198                                         console.log(sprintf('   \nLitho:'));
    199                                         fielddisplay(this,'numlayers','number of layers (default: 2)');
    200                                         fielddisplay(this,'radius','array describing the radius for each interface (numlayers+1) [m]');
    201                                         fielddisplay(this,'viscosity','array describing each layer\'s viscosity (numlayers) [Pa.s]');
    202                                         fielddisplay(this,'lame_lambda','array describing the lame lambda parameter (numlayers) [Pa]');
    203                                         fielddisplay(this,'lame_mu','array describing the shear modulus for each layers (numlayers) [Pa]');
    204                                         fielddisplay(this,'burgers_viscosity','array describing each layer\'s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]');
    205                                         fielddisplay(this,'burgers_mu','array describing each layer\'s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]');
    206 
    207                                         fielddisplay(this,'ebm_alpha','array describing each layer\'s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)');
    208                                         fielddisplay(this,'ebm_delta','array describing each layer\'s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)');
    209                                         fielddisplay(this,'ebm_taul','array describing each layer\'s starting period for transient relaxation, only for EBM rheology  (numlayers) [s]');
    210                                         fielddisplay(this,'ebm_tauh','array describing each layer\'s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]');
    211 
    212 
    213                                         fielddisplay(this,'rheologymodel','array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)');
    214                                         fielddisplay(this,'density','array describing each layer\'s density (numlayers) [kg/m^3]');
    215                                         fielddisplay(this,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)');
    216                                     break;
    217                                 case 'hydro':
    218                                         console.log(sprintf('   \nHydro:'));
    219                                         fielddisplay(this,'rho_ice','ice density [kg/m^3]');
    220                                         fielddisplay(this,'rho_water','ocean water density [kg/m^3]');
    221                                         fielddisplay(this,'earth_density','mantle density [kg/m^3]');
    222                                         fielddisplay(this,'rho_freshwater','fresh water density [kg/m^3]');
    223                     break;
    224                                 default:       
    225                                         error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
    226                         }       
    227                 }
    228 
    229         } // }}}
    230         extrude() {//{{{
    231                 for (let i=0; i < length(this.nature); i++) {
    232                         let nat=this.nature[i];
    233                         switch (nat) {
    234                                 case 'ice':
    235                                         this.rheology_B=project3d(md,'vector',this.rheology_B,'type','node');
    236                                         this.rheology_n=project3d(md,'vector',this.rheology_n,'type','element');
    237                                 break;
    238                         }
    239                 }
    240                 return this;
    241         }//}}}
    242         checkconsistency(md,solution,analyses) { // {{{
    243 
    244                 for (let i=0; i < length(this.nature); i++) {
    245                         let nat=this.nature[i];
    246                         switch (nat) {
    247                                 case 'ice':
    248                                         checkfield(md,'fieldname','materials.rho_ice','>',0);
    249                                         checkfield(md,'fieldname','materials.rho_water','>',0);
    250                                         checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    251                                         checkfield(md,'fieldname','materials.mu_water','>',0);
    252                                         checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1);
    253                                         checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements, 1]);
    254                                         checkfield(md,'fieldname','materials.rheology_law','values',['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O']);
    255                                     break;
    256                                 case 'litho':
    257                                         if (!ismember('LoveAnalysis',analyses)) return; end
    258                                         checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1);
    259                                         checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1, 1],'>',0);
    260                                         checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    261                                         checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    262                                         checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0,'<',2);
    263                                         checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>',0);
    264                                         checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    265                                         checkfield(md,'fieldname','materials.rheologymodel','NaN',1,'Inf',1,'size',[md.materials.numlayers, 1],'>=',0,'<=',2);
    266                                         checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    267                                         checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    268                                         checkfield(md,'fieldname','materials.ebm_alpha','Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    269                                         checkfield(md,'fieldname','materials.ebm_delta','Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    270                                         checkfield(md,'fieldname','materials.ebm_taul','Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    271                                         checkfield(md,'fieldname','materials.ebm_tauh','Inf',1,'size',[md.materials.numlayers, 1],'>=',0);
    272 
    273                                         for (let i=0; i < md.materials.numlayers; i++) {
    274                                                 if (md.materials.rheologymodel[i]==1 && (isNaN(md.materials.burgers_viscosity[i] || isNaN(md.materials.burgers_mu[i])))) {
    275                                                         error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice');
    276                                                 }
    277                                                 if (md.materials.rheologymodel[i]==2 && (isNaN(md.materials.ebm_alpha[i]) || isNaN(md.materials.ebm_delta[i]) || isNaN(md.materials.ebm_taul[i]) || isNaN(md.materials.ebm_tauh[i]))) {
    278                                                         error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice');
    279                                                 }
     166                }
     167        } //}}}
     168
     169        checkconsistency(md, solution, analyses) {//{{{
     170                for (let i = 0; i < this.nature.length; ++i) {
     171                        let nat = this.nature[i];
     172                        if (nat == 'ice') {
     173                                md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0);
     174                                md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0);
     175                                md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0);
     176                                md = checkfield(md, 'fieldname', 'materials.mu_water', '>', 0);
     177                                md = checkfield(md, 'fieldname', 'materials.rheology_B', '>', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1);
     178                                md = checkfield(md, 'fieldname', 'materials.rheology_n', '>', 0, 'size', [md.mesh.numberofelements, 1]);
     179                                md = checkfield(md, 'fieldname', 'materials.rheology_law', 'values', ['None', 'BuddJacka', 'Cuffey', 'CuffeyTemperate', 'Paterson', 'Arrhenius', 'LliboutryDuval', 'NyeCO2', 'NyeH2O']);
     180                        } else if (nat == 'litho') {
     181                                if (!analyses.includes('LoveAnalysis')) {
     182                                        return md;
     183                                }
     184
     185                                md = checkfield(md, 'fieldname', 'materials.numlayers', 'NaN', 1, 'Inf', 1, '>', 0, 'numel', 1);
     186                                md = checkfield(md, 'fieldname', 'materials.radius', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers + 1, 1], '>', 0);
     187                                md = checkfield(md, 'fieldname', 'materials.lame_mu', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     188                                md = checkfield(md, 'fieldname', 'materials.lame_lambda', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     189                                md = checkfield(md, 'fieldname', 'materials.issolid', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0, '<', 2);
     190                                md = checkfield(md, 'fieldname', 'materials.density', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>', 0);
     191                                md = checkfield(md, 'fieldname', 'materials.viscosity', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     192                                md = checkfield(md, 'fieldname', 'materials.rheologymodel', 'NaN', 1, 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0, '<=', 2);
     193                                md = checkfield(md, 'fieldname', 'materials.burgers_viscosity', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     194                                md = checkfield(md, 'fieldname', 'materials.burgers_mu', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     195                                md = checkfield(md, 'fieldname', 'materials.ebm_alpha', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     196                                md = checkfield(md, 'fieldname', 'materials.ebm_delta', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     197                                md = checkfield(md, 'fieldname', 'materials.ebm_taul', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     198                                md = checkfield(md, 'fieldname', 'materials.ebm_tauh', 'Inf', 1, 'size', [md.materials.numlayers, 1], '>=', 0);
     199
     200                                for (let i = 0; i < md.materials.numlayers; ++i) {
     201                                        if (md.materials.rheologymodel[i] == 1 && (isnan(md.materials.burgers_viscosity[i]) || isnan(md.materials.burgers_mu[i]))) {
     202                                                error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice');
    280203                                        }
    281                                         if (md.materials.issolid[0]==0 || md.materials.lame_mu[0]==0) {
    282                                                 error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.');
     204                                        if (md.materials.rheologymodel[i] == 2 && (isnan(md.materials.ebm_alpha[i]) || isnan(md.materials.ebm_delta[i])) || isnan(md.materials.ebm_taul[i] || isnan(md.materials.ebm_tauh[i]))) {
     205                                                error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice');
    283206                                        }
    284                                         ind=find(ArrayEqual(md.materials.issolid,0)); //ind=find(md.materials.issolid==0);
    285                                         console.log('materials consistency warning: not checking for fluid layers (not implemented)');
    286                                         //if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0
    287                                         //      error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
    288                                         //}
    289                     break;
    290                                 case 'hydro':
    291                                         checkfield(md,'fieldname','materials.rho_ice','>',0);
    292                                         checkfield(md,'fieldname','materials.rho_water','>',0);
    293                                         checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
    294                                         checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    295                     break;
    296                                 default:       
    297                                         error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
    298                         }
    299                 }
    300 
    301         } // }}}
    302         marshall(md,prefix,fid) { //{{{
    303                 //1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
    304                 WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(this.nature),'format','IntMat','mattype',3);
    305                 WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); //DANGER: this can evolve if you have classes.
    306                 for (let i=0; i < length(this.nature); i++) {
    307                         let nat=this.nature[i];
    308                         switch (nat) {
    309                                 case 'ice':
    310                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_ice','format','Double');
    311                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_water','format','Double');
    312                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_freshwater','format','Double');
    313                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','mu_water','format','Double');
    314                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','heatcapacity','format','Double');
    315                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','latentheat','format','Double');
    316                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','thermalconductivity','format','Double');
    317                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','temperateiceconductivity','format','Double');
    318                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','effectiveconductivity_averaging','format','Integer');
    319                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','meltingpoint','format','Double');
    320                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','beta','format','Double');
    321                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','mixed_layer_capacity','format','Double');
    322                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
    323                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    324                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
    325                                         WriteData(fid,prefix,'data',this.rheology_law,'name','md.materials.rheology_law','format','String');
    326                                     break;
    327                                 case 'litho':
    328                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','numlayers','format','Integer');
    329                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','radius','format','DoubleMat','mattype',3);
    330                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3);
    331                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3);
    332                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3);
    333                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','density','format','DoubleMat','mattype',3);
    334                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3);
    335                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rheologymodel','format','DoubleMat','mattype',3);
    336                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3);
    337                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3);
    338                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_alpha','format','DoubleMat','mattype',3);
    339                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_delta','format','DoubleMat','mattype',3);
    340                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_taul','format','DoubleMat','mattype',3);
    341                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','ebm_tauh','format','DoubleMat','mattype',3);
    342                                         //compute earth density compatible with our layer density distribution:
    343                                         let earth_density=0;
    344                                         for (let i=0; i < this.numlayers; i++){
    345                                                 earth_density=earth_density + (this.radius[i+1]**3-this.radius[i]**3)*this.density[i];
    346                                         }
    347                                         earth_density=earth_density/this.radius(this.numlayers+1)**3;
    348                                         console.log(earth_density);
    349                                         this.earth_density=earth_density;
    350                                     break;
    351                                 case 'hydro':
    352                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_ice','format','Double');
    353                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_water','format','Double');
    354                                         WriteData(fid,prefix,'object',this,'class','materials','fieldname','rho_freshwater','format','Double');
    355                                     break;
    356                                 default:       
    357                                         error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
    358                         }
    359                 }
    360                 WriteData(fid,prefix,'data',this.earth_density,'name','md.materials.earth_density','format','Double');
    361         }//}}}
    362 }
    363 function naturetointeger(strnat) {// {{{
    364         let intnat=zeros(length(strnat),1);
    365         for (let i=0; i < length(strnat); i++){
    366                 switch (strnat[i]) {
    367                         case 'damageice':
    368                                 intnat[i]=1;
    369                                 break;
    370                         case 'estar':
    371                                 intnat[i]=2;
    372                                 break;
    373                         case 'ice':
    374                                 intnat[i]=3;
    375                                 break;
    376                         case 'enhancedice':
    377                                 intnat[i]=4;
    378                                 break;
    379                         //case 'materials': //this case will never happen, kept to ensure equivalent of codes between IoCodeToMaterialsEnum and IoCodeToNatureEnum
    380                         //      intnat[i]=5;
    381                         //  break;
    382                         case 'litho':
    383                                 intnat[i]=6;
    384                                 break;
    385                         case 'hydro':
    386                                 intnat[i]=7;
    387                                 break;
    388                         default:
    389                                 error('materials constructor error message: nature of the material not supported yet! ("ice" or "litho" or "hydro")');
     207                                }
     208                                if (md.materials.issolid[0] == 0 || md.materials.lame_mu[0] == 0) {
     209                                        error('First layer must be solid (issolid[0] > 0 AND lame_mu[0] > 0). Add a weak inner core if necessary.');
     210                                }
     211                                let ind = find(md.materials.issolid == 0);
     212                                if (sum(ismember(diff(ind), 1) >= 1)) { // If there are at least two consecutive indices that contain issolid = 0
     213                                        error('Fluid layers detected at layers #' + ind + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.');
     214                                }
     215                        } else if (nat == 'hydro') {
     216                                md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0);
     217                                md = checkfield(md, 'fieldname', 'materials.rho_water', '>', 0);
     218                                md = checkfield(md, 'fieldname', 'materials.earth_density', '>', 0, 'numel', 1);
     219                                md = checkfield(md, 'fieldname', 'materials.rho_freshwater', '>', 0);
     220                        } else {
     221                                error('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')');
     222                        }
     223                }
     224
     225                return md;
     226        } //}}}
     227
     228        marshall(md, prefix, fid) {//{{{
     229                // 1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
     230                WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(this.nature), 'format', 'IntMat', 'mattype', 3);
     231                WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer'); // DANGER: this can evolve if you have classes
     232                for (let i = 0; i < this.nature.length; ++i) {
     233                        let nat = this.nature[i];
     234                        if (nat == 'ice') {
     235                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double');
     236                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double');
     237                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double');
     238                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'mu_water', 'format', 'Double');
     239                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'heatcapacity', 'format', 'Double');
     240                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'latentheat', 'format', 'Double');
     241                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'thermalconductivity', 'format', 'Double');
     242                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'temperateiceconductivity', 'format', 'Double');
     243                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'effectiveconductivity_averaging', 'format', 'Integer');
     244                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'meltingpoint', 'format', 'Double');
     245                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'beta', 'format', 'Double');
     246                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'mixed_layer_capacity', 'format', 'Double');
     247                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'thermal_exchange_velocity', 'format', 'Double');
     248                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rheology_B', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts);
     249                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rheology_n', 'format', 'DoubleMat', 'mattype', 2);
     250                                WriteData(fid, prefix, 'data', this.rheology_law, 'name', 'md.materials.rheology_law', 'format', 'String');
     251                        } else if (nat == 'litho') {
     252                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'numlayers', 'format', 'Integer');
     253                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'radius', 'format', 'DoubleMat', 'mattype', 3);
     254                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'lame_mu', 'format', 'DoubleMat', 'mattype', 3);
     255                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'lame_lambda', 'format', 'DoubleMat', 'mattype', 3);
     256                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'issolid', 'format', 'DoubleMat', 'mattype', 3);
     257                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'density', 'format', 'DoubleMat', 'mattype', 3);
     258                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'viscosity', 'format', 'DoubleMat', 'mattype', 3);
     259                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rheologymodel', 'format', 'DoubleMat', 'mattype', 3);
     260                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'burgers_viscosity', 'format', 'DoubleMat', 'mattype', 3);
     261                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'burgers_mu', 'format', 'DoubleMat', 'mattype', 3);
     262                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'ebm_alpha', 'format', 'DoubleMat', 'mattype', 3);
     263                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'ebm_delta', 'format', 'DoubleMat', 'mattype', 3);
     264                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'ebm_taul', 'format', 'DoubleMat', 'mattype', 3);
     265                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'ebm_tauh', 'format', 'DoubleMat', 'mattype', 3);
     266                                // Compute earth density compatible with our layer density distribution
     267                                let earth_density = 0;
     268                                for (let i = 0; i < this.numlayers; ++i) {
     269                                        earth_density = earth_density + (Math.pow(this.radius[i + 1], 3) - Math.pow(this.radius[i], 3)) * this.density[i];
     270                                }
     271                                earth_density = earth_density / Math.pow(this.radius[this.numlayers]);
     272                                this.earth_density = earth_density;
     273                        } else if (nat == 'hydro') {
     274                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double');
     275                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rho_water', 'format', 'Double');
     276                                WriteData(fid, prefix, 'object', this, 'class', 'materials', 'fieldname', 'rho_freshwater', 'format', 'Double');
     277                        } else {
     278                                error('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')');
     279                        }
     280                        WriteData(fid, prefix, 'data', this.earth_density, 'name', 'md.materials.earth_density', 'format', 'Double');
     281                }
     282        } //}}}
     283
     284        extrude(md) {//{{{
     285                for (let i = 0; i < this.nature.length; ++i) {
     286                        let nat = this.nature[i];
     287                        if (nat == 'ice') {
     288                                this.rheology_B = project3d(md, 'vector', this.rheology_B, 'type', 'node');
     289                                this.rheology_n = project3d(md, 'vector', this.rheology_n, 'type', 'element');
     290                        }
     291                }
     292        } //}}}
     293} //}}}
     294
     295function naturetointeger(strnat) {//{{{
     296        let intnat = zeros(strnat.length, 1);
     297        for (let i = 0; i < strnat.length; ++i) {
     298                let str_nat = strnat[i];
     299                if (str_nat == 'damageice') {
     300                        intnat[i] = 1;
     301                } else if (str_nat == 'estar') {
     302                        intnat[i] = 2;
     303                } else if (str_nat == 'ice') {
     304                        intnat[i] = 3;
     305                } else if (str_nat == 'enhancedice') {
     306                        intnat[i] = 4;
     307                //} else if (str_nat == 'materials') { // This case will never happen, kept to ensure equivalent of codes between IoCodeToMaterialsEnum and IoCodeToNatureEnum
     308                //      intnat[i] = 5;
     309                } else if (str_nat == 'litho') {
     310                        intnat[i] = 6;
     311                } else if (str_nat == 'hydro') {
     312                        intnat[i] = 7;
     313                } else {
     314                        error('materials constructor error message: nature of the material not supported yet! (\'ice\' or \'litho\' or \'hydro\')');
    390315                }
    391316        }
    392317        return intnat;
    393 }// }}}
     318} //}}}

  • issm/trunk-jpl/src/m/classes/materials.py

    r26301 r26308  
    332332
    333333    for i in range(len(intnat)):
    334         if strnat[i] == 'damageice':
     334        str_nat = strnat[i]
     335        if str_nat == 'damageice':
    335336            intnat[i] = 1
    336         elif strnat[i] == 'estar':
     337        elif str_nat == 'estar':
    337338            intnat[i] = 2
    338         elif strnat[i] == 'ice':
     339        elif str_nat == 'ice':
    339340            intnat[i] = 3
    340         elif strnat[i] == 'enhancedice':
     341        elif str_nat == 'enhancedice':
    341342            intnat[i] = 4
    342         # elif strnat[i] == 'materials':
     343        # elif str_nat == 'materials':
    343344        #     intnat[i] = 5 #this case will never happen, kept to ensure equivalent of codes between IoCodeToMeterialsEnum and IoCodeToNatureEnum
    344         elif strnat[i] == 'litho':
     345        elif str_nat == 'litho':
    345346            intnat[i] = 6
    346         elif strnat[i] == 'hydro':
     347        elif str_nat == 'hydro':
    347348            intnat[i] = 7
    348349        else:

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

    r26305 r26308  
    1 //MODEL class definition
    2 //
    3 //   Usage:
    4 //      md = new model()
    5 
    6 function model () {
     1/**
     2 * MODEL class definition
     3 *
     4 * Usage:
     5 *     md = new model();
     6 *
     7 * TODO:
     8 * - Convert to ES6 class
     9 */
     10function model(planet) {
    711        //methods
    8                 this.disp = function() { //{{{
    9                         console.log(sprintf("class model echo: "));
    10                         console.log(sprintf("//19s: //-22s -- //s","mesh"            ,"[1x1 " + typeof(this.mesh) + "]","mesh properties"));
    11                         console.log(sprintf("//19s: //-22s -- //s","mask"            ,"[1x1 " + typeof(this.mask) + "]","defines grounded and floating elements"));
    12                         console.log(sprintf("//19s: //-22s -- //s","geometry"        ,"[1x1 " + typeof(this.geometry) + "]","surface elevation, bedrock topography, ice thickness,..."));
    13                         console.log(sprintf("//19s: //-22s -- //s","constants"       ,"[1x1 " + typeof(this.constants) + "]","physical constants"));
    14                         console.log(sprintf("//19s: //-22s -- //s","smb"             ,"[1x1 " + typeof(this.smb) + "]","surface mass balance"));
    15                         console.log(sprintf("//19s: //-22s -- //s","basalforcings"   ,"[1x1 " + typeof(this.basalforcings) + "]","bed forcings"));
    16                         console.log(sprintf("//19s: //-22s -- //s","materials"       ,"[1x1 " + typeof(this.materials) + "]","material properties"));
    17                         console.log(sprintf("//19s: //-22s -- //s","damage"          ,"[1x1 " + typeof(this.damage) + "]","parameters for damage evolution solution"));
    18                         console.log(sprintf("//19s: //-22s -- //s","friction"        ,"[1x1 " + typeof(this.friction) + "]","basal friction/drag properties"));
    19                         console.log(sprintf("//19s: //-22s -- //s","flowequation"    ,"[1x1 " + typeof(this.flowequation) + "]","flow equations"));
    20                         console.log(sprintf("//19s: //-22s -- //s","timestepping"    ,"[1x1 " + typeof(this.timestepping) + "]","time stepping for trans models"));
    21                         console.log(sprintf("//19s: //-22s -- //s","initialization"  ,"[1x1 " + typeof(this.initialization) + "]","initial guess/state"));
    22                         console.log(sprintf("//19s: //-22s -- //s","rifts"           ,"[1x1 " + typeof(this.rifts) + "]","rifts properties"));
    23                         console.log(sprintf("//19s: //-22s -- //s","slr"             ,"[1x1 " + typeof(this.slr) + "]","slr forcings"));
    24                         console.log(sprintf("//19s: //-22s -- //s","dsl"             ,"[1x1 " + typeof(this.dsl) + "]","dynamic sea level"));
    25                         console.log(sprintf("//19s: //-22s -- //s","debug"           ,"[1x1 " + typeof(this.debug) + "]","debugging tools (valgrind, gprof)"));
    26                         console.log(sprintf("//19s: //-22s -- //s","verbose"         ,"[1x1 " + typeof(this.verbose) + "]","verbosity level in solve"));
    27                         console.log(sprintf("//19s: //-22s -- //s","settings"        ,"[1x1 " + typeof(this.settings) + "]","settings properties"));
    28                         console.log(sprintf("//19s: //-22s -- //s","toolkits"        ,"[1x1 " + typeof(this.toolkits) + "]","PETSc options for each solution"));
    29                         console.log(sprintf("//19s: //-22s -- //s","cluster"         ,"[1x1 " + typeof(this.cluster) + "]","cluster parameters (number of CPUs...)"));
    30                         console.log(sprintf("//19s: //-22s -- //s","balancethickness","[1x1 " + typeof(this.balancethickness) + "]","parameters for balancethickness solution"));
    31                         console.log(sprintf("//19s: //-22s -- //s","stressbalance"   ,"[1x1 " + typeof(this.stressbalance) + "]","parameters for stressbalance solution"));
    32                         console.log(sprintf("//19s: //-22s -- //s","groundingline"   ,"[1x1 " + typeof(this.groundingline) + "]","parameters for groundingline solution"));
    33                         console.log(sprintf("//19s: //-22s -- //s","hydrology"       ,"[1x1 " + typeof(this.hydrology) + "]","parameters for hydrology solution"));
    34                         console.log(sprintf("//19s: //-22s -- //s","masstransport"   ,"[1x1 " + typeof(this.masstransport) + "]","parameters for masstransport solution"));
    35                         console.log(sprintf("//19s: //-22s -- //s","thermal"         ,"[1x1 " + typeof(this.thermal) + "]","parameters for thermal solution"));
    36                         console.log(sprintf("//19s: //-22s -- //s","steadystate"     ,"[1x1 " + typeof(this.steadystate) + "]","parameters for steadystate solution"));
    37                         console.log(sprintf("//19s: //-22s -- //s","trans"           ,"[1x1 " + typeof(this.trans) + "]","parameters for trans solution"));
    38                         console.log(sprintf("//19s: //-22s -- //s","levelset"        ,"[1x1 " + typeof(this.levelset) + "]","parameters for moving boundaries (level-set method)"));
    39                         console.log(sprintf("//19s: //-22s -- //s","calving"         ,"[1x1 " + typeof(this.calving) + "]","parameters for calving"));
    40                         console.log(sprintf("//19s: //-22s -- //s","gia"             ,"[1x1 " + typeof(this.gia) + "]","parameters for gia solution"));
    41                         console.log(sprintf("//19s: //-22s -- //s","love"            ,"[1x1 " + typeof(this.love) + "]","parameters for love solution"));
    42                         console.log(sprintf("//19s: //-22s -- //s","esa"             ,"[1x1 " + typeof(this.esa) + "]","parameters for elastic adjustment solution"));
    43                         console.log(sprintf("//19s: //-22s -- //s","autodiff"        ,"[1x1 " + typeof(this.autodiff) + "]","automatic differentiation parameters"));
    44                         console.log(sprintf("//19s: //-22s -- //s","inversion"       ,"[1x1 " + typeof(this.inversion) + "]","parameters for inverse methods"));
    45                         console.log(sprintf("//19s: //-22s -- //s","qmu"             ,"[1x1 " + typeof(this.qmu) + "]","Dakota properties"));
    46                         console.log(sprintf("//19s: //-22s -- //s","amr"             ,"[1x1 " + typeof(this.amr) + "]","adaptive mesh refinement properties"));
    47                         console.log(sprintf("//19s: //-22s -- //s","outputdefinition","[1x1 " + typeof(this.outputdefinition) + "]","output definition"));
    48                         console.log(sprintf("//19s: //-22s -- //s","results"         ,"[1x1 " + typeof(this.results) + "]","model results"));
    49                         console.log(sprintf("//19s: //-22s -- //s","radaroverlay"    ,"[1x1 " + typeof(this.radaroverlay) + "]","radar image for plot overlay"));
    50                         console.log(sprintf("//19s: //-22s -- //s","miscellaneous"   ,"[1x1 " + typeof(this.miscellaneous) + "]","miscellaneous fields"));
    51                 } //}}}
    52                 this.setdefaultparameters = function () { // {{{
    53 
    54                         //initialize subclasses
    55                         this.mesh             = new mesh2d();
    56                         this.mask             = new mask();
    57                         this.constants        = new constants();
    58                         this.geometry         = new geometry();
    59                         this.smb              = new SMBforcing();
    60                         this.basalforcings    = new basalforcings();
    61                         this.materials        = new matice();
    62                         this.damage           = new damage();
    63                         this.friction         = new friction();
    64                         this.flowequation     = new flowequation();
    65                         this.timestepping     = new timestepping();
    66                         this.initialization   = new initialization();
    67                         this.rifts            = new rifts();
    68                         this.slr              = new slr();
    69                         this.dsl              = new dsl();
    70 
    71                         this.debug            = new debug();
    72                         this.verbose          = new verbose();
    73                         this.settings         = new issmsettings();
    74                         this.toolkits         = new toolkits();
    75                         this.cluster          = new local();
    76 
    77                         this.balancethickness = new balancethickness();
    78                         this.stressbalance    = new stressbalance();
    79                         this.groundingline    = new groundingline();
    80                         this.hydrology        = new hydrologyshreve();
    81                         this.masstransport    = new masstransport();
    82                         this.thermal          = new thermal();
    83                         this.steadystate      = new steadystate();
    84                         this.trans            = new trans();
    85                         this.levelset         = new levelset();
    86                         this.calving          = new calving();
    87                         this.gia              = new giaivins();
    88                         this.frontforcings    = new frontforcings();
    89                         this.love             = new fourierlove();
    90                         this.esa              = new esa();
    91                         this.autodiff         = new autodiff();
    92                         this.inversion        = new inversion();
    93                         this.qmu              = new qmu();
    94                         this.amr                                         = new amr();
    95                         this.radaroverlay     = new radaroverlay();
    96                         this.results          = {};
    97                         this.outputdefinition = new outputdefinition();
    98                         this.radaroverlay     = new radaroverlay();
    99                         this.miscellaneous    = new miscellaneous();
    100                         this.priv             = new priv();
    101                 } //}}}
    102                 this.checkmessage = function(string){ //{{{
    103                         console.log('model not consistent: ' + string);
    104                         this.priv.isconsistent=false;
    105                 } //}}}
    106                 this.fix = function(){ //{{{
    107 
    108                         for (var field in this){
    109 
    110                                 //Some properties do not need to be fixed
    111                                 if (field == 'results' | field =='radaroverlay' | field == 'toolkits' | field =='cluster' | field == 'priv') continue;
    112 
    113                                 //Check that current field is a class
    114                                 if(typeof this[field] == 'function'){
    115                                         continue;
     12        this.disp = function() { //{{{
     13                console.log(sprintf("class model echo: "));
     14                console.log(sprintf("//19s: //-22s -- //s","mesh"            ,"[1x1 " + typeof(this.mesh) + "]","mesh properties"));
     15                console.log(sprintf("//19s: //-22s -- //s","mask"            ,"[1x1 " + typeof(this.mask) + "]","defines grounded and floating elements"));
     16                console.log(sprintf("//19s: //-22s -- //s","geometry"        ,"[1x1 " + typeof(this.geometry) + "]","surface elevation, bedrock topography, ice thickness,..."));
     17                console.log(sprintf("//19s: //-22s -- //s","constants"       ,"[1x1 " + typeof(this.constants) + "]","physical constants"));
     18                console.log(sprintf("//19s: //-22s -- //s","smb"             ,"[1x1 " + typeof(this.smb) + "]","surface mass balance"));
     19                console.log(sprintf("//19s: //-22s -- //s","basalforcings"   ,"[1x1 " + typeof(this.basalforcings) + "]","bed forcings"));
     20                console.log(sprintf("//19s: //-22s -- //s","materials"       ,"[1x1 " + typeof(this.materials) + "]","material properties"));
     21                console.log(sprintf("//19s: //-22s -- //s","damage"          ,"[1x1 " + typeof(this.damage) + "]","parameters for damage evolution solution"));
     22                console.log(sprintf("//19s: //-22s -- //s","friction"        ,"[1x1 " + typeof(this.friction) + "]","basal friction/drag properties"));
     23                console.log(sprintf("//19s: //-22s -- //s","flowequation"    ,"[1x1 " + typeof(this.flowequation) + "]","flow equations"));
     24                console.log(sprintf("//19s: //-22s -- //s","timestepping"    ,"[1x1 " + typeof(this.timestepping) + "]","time stepping for trans models"));
     25                console.log(sprintf("//19s: //-22s -- //s","initialization"  ,"[1x1 " + typeof(this.initialization) + "]","initial guess/state"));
     26                console.log(sprintf("//19s: //-22s -- //s","rifts"           ,"[1x1 " + typeof(this.rifts) + "]","rifts properties"));
     27                console.log(sprintf("//19s: //-22s -- //s","solidearth"      ,"[1x1 " + typeof(this.solidearth) + "]","solidearth inputs and settings"));
     28                console.log(sprintf("//19s: //-22s -- //s","slr"             ,"[1x1 " + typeof(this.slr) + "]","slr forcings"));
     29                console.log(sprintf("//19s: //-22s -- //s","dsl"             ,"[1x1 " + typeof(this.dsl) + "]","dynamic sea level"));
     30                console.log(sprintf("//19s: //-22s -- //s","debug"           ,"[1x1 " + typeof(this.debug) + "]","debugging tools (valgrind, gprof)"));
     31                console.log(sprintf("//19s: //-22s -- //s","verbose"         ,"[1x1 " + typeof(this.verbose) + "]","verbosity level in solve"));
     32                console.log(sprintf("//19s: //-22s -- //s","settings"        ,"[1x1 " + typeof(this.settings) + "]","settings properties"));
     33                console.log(sprintf("//19s: //-22s -- //s","toolkits"        ,"[1x1 " + typeof(this.toolkits) + "]","PETSc options for each solution"));
     34                console.log(sprintf("//19s: //-22s -- //s","cluster"         ,"[1x1 " + typeof(this.cluster) + "]","cluster parameters (number of CPUs...)"));
     35                console.log(sprintf("//19s: //-22s -- //s","balancethickness","[1x1 " + typeof(this.balancethickness) + "]","parameters for balancethickness solution"));
     36                console.log(sprintf("//19s: //-22s -- //s","stressbalance"   ,"[1x1 " + typeof(this.stressbalance) + "]","parameters for stressbalance solution"));
     37                console.log(sprintf("//19s: //-22s -- //s","groundingline"   ,"[1x1 " + typeof(this.groundingline) + "]","parameters for groundingline solution"));
     38                console.log(sprintf("//19s: //-22s -- //s","hydrology"       ,"[1x1 " + typeof(this.hydrology) + "]","parameters for hydrology solution"));
     39                console.log(sprintf("//19s: //-22s -- //s","masstransport"   ,"[1x1 " + typeof(this.masstransport) + "]","parameters for masstransport solution"));
     40                console.log(sprintf("//19s: //-22s -- //s","thermal"         ,"[1x1 " + typeof(this.thermal) + "]","parameters for thermal solution"));
     41                console.log(sprintf("//19s: //-22s -- //s","steadystate"     ,"[1x1 " + typeof(this.steadystate) + "]","parameters for steadystate solution"));
     42                console.log(sprintf("//19s: //-22s -- //s","transient"       ,"[1x1 " + typeof(this.transient) + "]","parameters for transient solution"));
     43                console.log(sprintf("//19s: //-22s -- //s","levelset"        ,"[1x1 " + typeof(this.levelset) + "]","parameters for moving boundaries (level-set method)"));
     44                console.log(sprintf("//19s: //-22s -- //s","calving"         ,"[1x1 " + typeof(this.calving) + "]","parameters for calving"));
     45                console.log(sprintf("//19s: //-22s -- //s","gia"             ,"[1x1 " + typeof(this.gia) + "]","parameters for gia solution"));
     46                console.log(sprintf("//19s: //-22s -- //s","love"            ,"[1x1 " + typeof(this.love) + "]","parameters for love solution"));
     47                console.log(sprintf("//19s: //-22s -- //s","esa"             ,"[1x1 " + typeof(this.esa) + "]","parameters for elastic adjustment solution"));
     48                console.log(sprintf("//19s: //-22s -- //s","sampling"        ,"[1x1 " + typeof(this.sampling) + "]","parameters for stochastic sampler"));
     49                console.log(sprintf("//19s: //-22s -- //s","autodiff"        ,"[1x1 " + typeof(this.autodiff) + "]","automatic differentiation parameters"));
     50                console.log(sprintf("//19s: //-22s -- //s","inversion"       ,"[1x1 " + typeof(this.inversion) + "]","parameters for inverse methods"));
     51                console.log(sprintf("//19s: //-22s -- //s","qmu"             ,"[1x1 " + typeof(this.qmu) + "]","Dakota properties"));
     52                console.log(sprintf("//19s: //-22s -- //s","amr"             ,"[1x1 " + typeof(this.amr) + "]","adaptive mesh refinement properties"));
     53                console.log(sprintf("//19s: //-22s -- //s","outputdefinition","[1x1 " + typeof(this.outputdefinition) + "]","output definition"));
     54                console.log(sprintf("//19s: //-22s -- //s","results"         ,"[1x1 " + typeof(this.results) + "]","model results"));
     55                console.log(sprintf("//19s: //-22s -- //s","radaroverlay"    ,"[1x1 " + typeof(this.radaroverlay) + "]","radar image for plot overlay"));
     56                console.log(sprintf("//19s: //-22s -- //s","miscellaneous"   ,"[1x1 " + typeof(this.miscellaneous) + "]","miscellaneous fields"));
     57        } //}}}
     58        this.setdefaultparameters = function(planet) { // {{{
     59                this.mesh             = new mesh2d();
     60                this.mask             = new mask();
     61
     62                this.geometry         = new geometry();
     63                this.constants        = new constants();
     64                this.smb              = new SMBforcing();
     65                this.basalforcings    = new basalforcings();
     66                this.materials        = new matice();
     67                this.damage           = new damage();
     68                this.friction         = new friction();
     69                this.flowequation     = new flowequation();
     70                this.timestepping     = new timestepping();
     71                this.initialization   = new initialization();
     72                this.rifts            = new rifts();
     73                this.dsl              = new dsl();
     74                this.solidearth       = new solidearth(planet);
     75
     76                this.debug            = new debug();
     77                this.verbose          = new verbose();
     78                this.settings         = new issmsettings();
     79                this.toolkits         = new toolkits();
     80                this.cluster          = new local();
     81
     82                this.balancethickness = new balancethickness();
     83                this.stressbalance    = new stressbalance();
     84                this.groundingline    = new groundingline();
     85                this.hydrology        = new hydrologyshreve();
     86                this.masstransport    = new masstransport();
     87                this.thermal          = new thermal();
     88                this.steadystate      = new steadystate();
     89                this.transient        = new transient();
     90                this.levelset         = new levelset();
     91                this.calving          = new calving();
     92                this.frontalforcings  = new frontalforcings();
     93                this.love             = new fourierlove();
     94                this.esa              = new esa();
     95                this.sampling         = new sampling();
     96
     97                this.autodiff         = new autodiff();
     98                this.inversion        = new inversion();
     99                this.qmu              = new qmu();
     100                this.amr              = new amr();
     101                this.results          = {};
     102                this.outputdefinition = new outputdefinition();
     103                this.radaroverlay     = new radaroverlay();
     104                this.miscellaneous    = new miscellaneous();
     105                this.priv             = new priv();
     106        } //}}}
     107        this.checkmessage = function(string){ //{{{
     108                console.log('model not consistent: ' + string);
     109                this.priv.isconsistent=false;
     110        } //}}}
     111        this.fix = function(){ //{{{
     112
     113                for (var field in this){
     114
     115                        //Some properties do not need to be fixed
     116                        if (field == 'results' | field =='radaroverlay' | field == 'toolkits' | field =='cluster' | field == 'priv') continue;
     117
     118                        //Check that current field is a class
     119                        if(typeof this[field] == 'function'){
     120                                continue;
     121                        }
     122
     123                        //Fix current object
     124                        this[field].fix(this);
     125                }
     126
     127        } //}}}
     128        this.extrude = function(md) { //{{{
     129                //EXTRUDE - vertically extrude a 2d mesh
     130                //
     131                //   vertically extrude a 2d mesh and create corresponding 3d mesh.
     132                //   The vertical distribution can:
     133                //    - follow a polynomial law
     134                //    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
     135                //    - be described by a list of coefficients (between 0 and 1)
     136                //   
     137                //
     138                //   Usage:
     139                //      md=extrude(md,numlayers,extrusionexponent);
     140                //      md=extrude(md,numlayers,lowerexponent,upperexponent);
     141                //      md=extrude(md,listofcoefficients);
     142                //
     143                //   Example:
     144                //      md=extrude(md,15,1.3);
     145                //      md=extrude(md,15,1.3,1.2);
     146                //      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
     147                //
     148                //   See also: MODELEXTRACT, COLLAPSE
     149
     150                //some checks on list of arguments
     151                var argc = arguments.length;
     152                var extrusionlist;
     153                var numlayers;
     154               
     155                if ((argc > 4) | (argc < 2)) {
     156                        error("extrude error message");
     157                }
     158
     159                //Extrude the mesh
     160                if (argc==2) { //list of coefficients
     161                        clist=arguments[0];
     162                        if (ArrayAnyBelowStrict(clist,0) || ArrayAnyAboveStrict(clist,1)) {
     163                                error('extrusioncoefficients must be between 0 and 1');
     164                        }
     165                        extrusionlist=ArraySort(ArrayUnique(clist.push(0,1)));
     166                        numlayers=extrusionlist.length;
     167                } else if (argc==3) { //one polynomial law
     168                        if (arguments[1]<=0) {
     169                                error('extrusionexponent must be >=0');
     170                        }
     171                        numlayers=arguments[1];
     172                        extrusionlist = [];
     173                        for (var i = 0; i < numlayers; i++) {
     174                                extrusionlist.push(Math.pow(i/(numlayers-1), arguments[2]));
     175                        }
     176                } else if (argc==4) { //two polynomial laws
     177                        numlayers=arguments[1];
     178                        var lowerexp=arguments[2];
     179                        var upperexp=arguments[3];
     180
     181                        if (arguments[2]<=0 || arguments[3]<=0) {
     182                                error('lower and upper extrusionexponents must be >=0');
     183                        }
     184
     185                        var lowerextrusionlist = [];
     186                        for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
     187                                lowerextrusionlist.push(Math.pow(i, lowerexp)/2)
     188                        }
     189                        var upperextrusionlist = [];
     190                        for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
     191                                upperextrusionlist.push(1-Math.pow(i, upperexp)/2)
     192                        }
     193                        extrusionlist=ArrayConcat(lowerextrusionlist,upperextrusionlist);
     194                        extrusionlist=ArraySort(ArrayUnique(extrusionlist.push(1)));
     195                }
     196
     197                if (numlayers<2) {
     198                        console.error('number of layers should be at least 2');
     199                }
     200                if (md.mesh.domaintype() === '3D') {
     201                        console.error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
     202                }
     203
     204                //Initialize with the 2d mesh
     205                var mesh2d = md.mesh;
     206                md.mesh=new mesh3dprisms();
     207                md.mesh.x                           = mesh2d.x;
     208                md.mesh.y                           = mesh2d.y;
     209                md.mesh.elements                    = mesh2d.elements;
     210                md.mesh.numberofelements            = mesh2d.numberofelements;
     211                md.mesh.numberofvertices            = mesh2d.numberofvertices;
     212
     213                md.mesh.lat                         = mesh2d.lat;
     214                md.mesh.long                        = mesh2d.long;
     215                md.mesh.epsg                        = mesh2d.epsg;
     216                md.mesh.scale_factor                = mesh2d.scale_factor;
     217
     218                md.mesh.vertexonboundary            = mesh2d.vertexonboundary;
     219                md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity;
     220                md.mesh.elementconnectivity         = mesh2d.elementconnectivity;
     221                md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
     222
     223                md.mesh.extractedvertices           = mesh2d.extractedvertices;
     224                md.mesh.extractedelements           = mesh2d.extractedelements;
     225
     226                var x3d=new Float64Array();
     227                var y3d=new Float64Array();
     228                var z3d=new Float64Array();  //the lower node is on the bed
     229                var thickness3d=md.geometry.thickness; //thickness and bed for these nodes
     230                var bed3d=md.geometry.base;
     231
     232                //Create the new layers
     233                //Dimensions of x/y/z3d: md.mesh.numberofvertices * numlayers, 1
     234                for (var i = 1; i <= numlayers; i++) {
     235                        x3d=ArrayConcat(x3d, md.mesh.x);
     236                        y3d=ArrayConcat(y3d, md.mesh.y);
     237                        //nodes are distributed between bed and surface accordingly to the given exponent
     238                        z3d=ArrayConcat(z3d, FloatFix(thickness3d.map(function(value,index) { return bed3d[index] + value * extrusionlist[i-1]; }), thickness3d.length));
     239                }
     240                var number_nodes3d=x3d.length; //number of 3d nodes for the non extruded part of the mesh
     241
     242                //Extrude elements
     243                //Create the elements of the 3d mesh for the non extruded part
     244                //Dimensions of elements3d: md.mesh.numberofelements * (numlayers - 1), 6
     245                var elements3d=[];
     246                var elements_prisms=[];
     247                for (var i = 1; i < numlayers; i++) {
     248                        for (var j = 0; j < md.mesh.numberofelements; j++) {
     249                                elements_prisms = [];
     250                                for (var k = 0, offset = (i - 1) * md.mesh.numberofvertices; k < 3; k++) {
     251                                        elements_prisms.push(md.mesh.elements[j][k]+offset);
    116252                                }
    117 
    118                                 //Fix current object
    119                                 this[field].fix(this);
    120                         }
    121 
    122                 } //}}}
    123                 this.extrude = function(md) { //{{{
    124                         //EXTRUDE - vertically extrude a 2d mesh
    125                         //
    126                         //   vertically extrude a 2d mesh and create corresponding 3d mesh.
    127                         //   The vertical distribution can:
    128                         //    - follow a polynomial law
    129                         //    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
    130                         //    - be described by a list of coefficients (between 0 and 1)
    131                         //   
    132                         //
    133                         //   Usage:
    134                         //      md=extrude(md,numlayers,extrusionexponent);
    135                         //      md=extrude(md,numlayers,lowerexponent,upperexponent);
    136                         //      md=extrude(md,listofcoefficients);
    137                         //
    138                         //   Example:
    139                         //      md=extrude(md,15,1.3);
    140                         //      md=extrude(md,15,1.3,1.2);
    141                         //      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
    142                         //
    143                         //   See also: MODELEXTRACT, COLLAPSE
    144 
    145                         //some checks on list of arguments
    146                         var argc = arguments.length;
    147                         var extrusionlist;
    148                         var numlayers;
    149                        
    150                         if ((argc > 4) | (argc < 2)) {
    151                                 throw "extrude error message";
    152                         }
    153 
    154                         //Extrude the mesh
    155                         if (argc==2) { //list of coefficients
    156                                 clist=arguments[0];
    157                                 if (ArrayAnyBelowStrict(clist,0) || ArrayAnyAboveStrict(clist,1)) {
    158                                         throw 'extrusioncoefficients must be between 0 and 1';
     253                                for (var k = 0, offset = i * md.mesh.numberofvertices; k < 3; k++) {
     254                                        elements_prisms.push(md.mesh.elements[j][k]+offset);
    159255                                }
    160                                 extrusionlist=ArraySort(ArrayUnique(clist.push(0,1)));
    161                                 numlayers=extrusionlist.length;
    162                         } else if (argc==3) { //one polynomial law
    163                                 if (arguments[1]<=0) {
    164                                         throw 'extrusionexponent must be >=0';
     256                                elements3d.push(elements_prisms);
     257                        }                               
     258                }
     259                number_el3d=elements3d.length; //number of 3d nodes for the non extruded part of the mesh
     260
     261                //Keep a trace of lower and upper nodes
     262                var lowervertex = NewArrayFill(number_nodes3d, NaN);
     263                var uppervertex = NewArrayFill(number_nodes3d, NaN);
     264                //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofvertices+1)
     265                for (var i = md.mesh.numberofvertices, k = 1; i < lowervertex.length; i++, k++) {
     266                        lowervertex[i] = k;
     267                };
     268                //Set last layer to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1)
     269                for (var i = 0, k = md.mesh.numberofvertices+1; i < (numlayers-1)*md.mesh.numberofvertices; i++, k++) {
     270                        uppervertex[i] = k;
     271                };
     272                md.mesh.lowervertex=lowervertex;
     273                md.mesh.uppervertex=uppervertex;
     274
     275                //same for lower and upper elements
     276                var lowerelements = NewArrayFill(number_el3d, NaN);
     277                var upperelements = NewArrayFill(number_el3d, NaN);
     278                //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofelements+1)
     279                for (var i = md.mesh.numberofelements, k = 1; i < lowerelements.length; i++, k++) {
     280                        lowerelements[i] = k;
     281                };
     282                //Set last 2 layers to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1)
     283                for (var i = 0, k = md.mesh.numberofelements + 1; i < (numlayers-2)*md.mesh.numberofelements; i++, k++) {
     284                        upperelements[i] = k;
     285                };
     286                md.mesh.lowerelements=lowerelements;
     287                md.mesh.upperelements=upperelements;
     288
     289                //Save old mesh
     290                md.mesh.x2d=md.mesh.x;
     291                md.mesh.y2d=md.mesh.y;
     292                md.mesh.elements2d=md.mesh.elements;
     293                md.mesh.numberofelements2d=md.mesh.numberofelements;
     294                md.mesh.numberofvertices2d=md.mesh.numberofvertices;
     295
     296                //Build global 3d mesh
     297                md.mesh.elements=elements3d;
     298                md.mesh.x=x3d;
     299                md.mesh.y=y3d;
     300                md.mesh.z=z3d;
     301                md.mesh.numberofelements=number_el3d;
     302                md.mesh.numberofvertices=number_nodes3d;
     303                md.mesh.numberoflayers=numlayers;
     304
     305                //Ok, now deal with the other fields from the 2d mesh:
     306
     307                //bedinfo and surface info
     308                md.mesh.vertexonbase=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',1);
     309                md.mesh.vertexonsurface=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
     310                md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
     311
     312                //lat long
     313                md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
     314                md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
     315                md.mesh.scale_factor=project3d(md,'vector',md.mesh.scale_factor,'type','node');
     316
     317                md.geometry=md.geometry.extrude(md);
     318                md.friction=md.friction.extrude(md);
     319                md.inversion=md.inversion.extrude(md);
     320                md.smb=md.smb.extrude(md);
     321                md.initialization=md.initialization.extrude(md);
     322
     323                md.flowequation=md.flowequation.extrude(md);
     324                md.stressbalance=md.stressbalance.extrude(md);
     325                md.thermal=md.thermal.extrude(md);
     326                md.masstransport=md.masstransport.extrude(md);
     327                md.levelset=md.levelset.extrude(md);
     328                md.calving=md.calving.extrude(md);
     329                md.hydrology = md.hydrology.extrude(md);
     330                md.solidearth = md.solidearth.extrude(md);
     331
     332                //connectivity
     333                if (!ArrayAnyNaN(md.mesh.elementconnectivity)) {
     334                        //md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); //Replicate the matrix across numlayers-1
     335                        var temparr = [];
     336                        for (var i = 0; i < numlayers - 1; i++) {
     337                                temparr = ArrayConcat(temparr, md.mesh.elementconnectivity);
     338                        }
     339                        md.mesh.elementconnectivity = temparr;
     340       
     341                        //md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
     342                        var indices = ArrayFind(md.mesh.elementconnectivity, 0);
     343                        for (var i = 0; i < indices.length; i++) {
     344                                md.mesh.elementconnectivity[i] = NaN;
     345                        };
     346       
     347                        for (var i = 2; i < numlayers; i++) {
     348                                //md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
     349                                //=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
     350                                for (var j = (i-1)*md.mesh.numberofelements2d; j <= (i)*md.mesh.numberofelements2d-1; j++) {
     351                                        for (var k = 0; k < 3; k++) {
     352                                                md.mesh.elementconnectivity[j][k] += md.mesh.numberofelements2d;
     353                                        }
    165354                                }
    166                                 numlayers=arguments[1];
    167                                 extrusionlist = [];
    168                                 for (var i = 0; i < numlayers; i++) {
    169                                         extrusionlist.push(Math.pow(i/(numlayers-1), arguments[2]));
    170                                 }
    171                         } else if (argc==4) { //two polynomial laws
    172                                 numlayers=arguments[1];
    173                                 var lowerexp=arguments[2];
    174                                 var upperexp=arguments[3];
    175 
    176                                 if (arguments[2]<=0 || arguments[3]<=0) {
    177                                         throw 'lower and upper extrusionexponents must be >=0';
    178                                 }
    179 
    180                                 var lowerextrusionlist = [];
    181                                 for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
    182                                         lowerextrusionlist.push(Math.pow(i, lowerexp)/2)
    183                                 }
    184                                 var upperextrusionlist = [];
    185                                 for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
    186                                         upperextrusionlist.push(1-Math.pow(i, upperexp)/2)
    187                                 }
    188                                 extrusionlist=ArrayConcat(lowerextrusionlist,upperextrusionlist);
    189                                 extrusionlist=ArraySort(ArrayUnique(extrusionlist.push(1)));
    190                         }
    191 
    192                         if (numlayers<2) {
    193                                 console.error('number of layers should be at least 2');
    194                         }
    195                         if (md.mesh.domaintype() === '3D') {
    196                                 console.error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
    197                         }
    198 
    199                         //Initialize with the 2d mesh
    200                         var mesh2d = md.mesh;
    201                         md.mesh=new mesh3dprisms();
    202                         md.mesh.x                           = mesh2d.x;
    203                         md.mesh.y                           = mesh2d.y;
    204                         md.mesh.elements                    = mesh2d.elements;
    205                         md.mesh.numberofelements            = mesh2d.numberofelements;
    206                         md.mesh.numberofvertices            = mesh2d.numberofvertices;
    207 
    208                         md.mesh.lat                         = mesh2d.lat;
    209                         md.mesh.long                        = mesh2d.long;
    210                         md.mesh.epsg                        = mesh2d.epsg;
    211                         md.mesh.scale_factor                = mesh2d.scale_factor;
    212 
    213                         md.mesh.vertexonboundary            = mesh2d.vertexonboundary;
    214                         md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity;
    215                         md.mesh.elementconnectivity         = mesh2d.elementconnectivity;
    216                         md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;
    217 
    218                         md.mesh.extractedvertices           = mesh2d.extractedvertices;
    219                         md.mesh.extractedelements           = mesh2d.extractedelements;
    220 
    221                         var x3d=new Float64Array();
    222                         var y3d=new Float64Array();
    223                         var z3d=new Float64Array();  //the lower node is on the bed
    224                         var thickness3d=md.geometry.thickness; //thickness and bed for these nodes
    225                         var bed3d=md.geometry.base;
    226 
    227                         //Create the new layers
    228                         //Dimensions of x/y/z3d: md.mesh.numberofvertices * numlayers, 1
    229                         for (var i = 1; i <= numlayers; i++) {
    230                                 x3d=ArrayConcat(x3d, md.mesh.x);
    231                                 y3d=ArrayConcat(y3d, md.mesh.y);
    232                                 //nodes are distributed between bed and surface accordingly to the given exponent
    233                                 z3d=ArrayConcat(z3d, FloatFix(thickness3d.map(function(value,index) { return bed3d[index] + value * extrusionlist[i-1]; }), thickness3d.length));
    234                         }
    235                         var number_nodes3d=x3d.length; //number of 3d nodes for the non extruded part of the mesh
    236 
    237                         //Extrude elements
    238                         //Create the elements of the 3d mesh for the non extruded part
    239                         //Dimensions of elements3d: md.mesh.numberofelements * (numlayers - 1), 6
    240                         var elements3d=[];
    241                         var elements_prisms=[];
    242                         for (var i = 1; i < numlayers; i++) {
    243                                 for (var j = 0; j < md.mesh.numberofelements; j++) {
    244                                         elements_prisms = [];
    245                                         for (var k = 0, offset = (i - 1) * md.mesh.numberofvertices; k < 3; k++) {
    246                                                 elements_prisms.push(md.mesh.elements[j][k]+offset);
    247                                         }
    248                                         for (var k = 0, offset = i * md.mesh.numberofvertices; k < 3; k++) {
    249                                                 elements_prisms.push(md.mesh.elements[j][k]+offset);
    250                                         }
    251                                         elements3d.push(elements_prisms);
    252                                 }                               
    253                         }
    254                         number_el3d=elements3d.length; //number of 3d nodes for the non extruded part of the mesh
    255 
    256                         //Keep a trace of lower and upper nodes
    257                         var lowervertex = NewArrayFill(number_nodes3d, NaN);
    258                         var uppervertex = NewArrayFill(number_nodes3d, NaN);
    259                         //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofvertices+1)
    260                         for (var i = md.mesh.numberofvertices, k = 1; i < lowervertex.length; i++, k++) {
    261                                 lowervertex[i] = k;
    262                         };
    263                         //Set last layer to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1)
    264                         for (var i = 0, k = md.mesh.numberofvertices+1; i < (numlayers-1)*md.mesh.numberofvertices; i++, k++) {
    265                                 uppervertex[i] = k;
    266                         };
    267                         md.mesh.lowervertex=lowervertex;
    268                         md.mesh.uppervertex=uppervertex;
    269 
    270                         //same for lower and upper elements
    271                         var lowerelements = NewArrayFill(number_el3d, NaN);
    272                         var upperelements = NewArrayFill(number_el3d, NaN);
    273                         //Set first layer to NaN and start count from 1 at next layer (i = md.mesh.numberofelements+1)
    274                         for (var i = md.mesh.numberofelements, k = 1; i < lowerelements.length; i++, k++) {
    275                                 lowerelements[i] = k;
    276                         };
    277                         //Set last 2 layers to NaN and start count from md.mesh.numberofvertices+1 at first layer (i = 1)
    278                         for (var i = 0, k = md.mesh.numberofelements + 1; i < (numlayers-2)*md.mesh.numberofelements; i++, k++) {
    279                                 upperelements[i] = k;
    280                         };
    281                         md.mesh.lowerelements=lowerelements;
    282                         md.mesh.upperelements=upperelements;
    283 
    284                         //Save old mesh
    285                         md.mesh.x2d=md.mesh.x;
    286                         md.mesh.y2d=md.mesh.y;
    287                         md.mesh.elements2d=md.mesh.elements;
    288                         md.mesh.numberofelements2d=md.mesh.numberofelements;
    289                         md.mesh.numberofvertices2d=md.mesh.numberofvertices;
    290 
    291                         //Build global 3d mesh
    292                         md.mesh.elements=elements3d;
    293                         md.mesh.x=x3d;
    294                         md.mesh.y=y3d;
    295                         md.mesh.z=z3d;
    296                         md.mesh.numberofelements=number_el3d;
    297                         md.mesh.numberofvertices=number_nodes3d;
    298                         md.mesh.numberoflayers=numlayers;
    299 
    300                         //Ok, now deal with the other fields from the 2d mesh:
    301 
    302                         //bedinfo and surface info
    303                         md.mesh.vertexonbase=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',1);
    304                         md.mesh.vertexonsurface=project3d(md,'vector',NewArrayFill(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
    305                         md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');
    306 
    307                         //lat long
    308                         md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
    309                         md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
    310                         md.mesh.scale_factor=project3d(md,'vector',md.mesh.scale_factor,'type','node');
    311 
    312                         md.geometry=md.geometry.extrude(md);
    313                         md.friction=md.friction.extrude(md);
    314                         md.inversion=md.inversion.extrude(md);
    315                         md.smb=md.smb.extrude(md);
    316                         md.initialization=md.initialization.extrude(md);
    317 
    318                         md.flowequation=md.flowequation.extrude(md);
    319                         md.stressbalance=md.stressbalance.extrude(md);
    320                         md.thermal=md.thermal.extrude(md);
    321                         md.masstransport=md.masstransport.extrude(md);
    322                         md.levelset=md.levelset.extrude(md);
    323                         md.calving=md.calving.extrude(md);
    324                         md.hydrology = md.hydrology.extrude(md);
    325 
    326                         //connectivity
    327                         if (!ArrayAnyNaN(md.mesh.elementconnectivity)) {
    328                                 //md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1); //Replicate the matrix across numlayers-1
    329                                 var temparr = [];
    330                                 for (var i = 0; i < numlayers - 1; i++) {
    331                                         temparr = ArrayConcat(temparr, md.mesh.elementconnectivity);
    332                                 }
    333                                 md.mesh.elementconnectivity = temparr;
     355                        }
     356                        md.mesh.elementconnectivity = md.mesh.elementconnectivity.map(function(value) { return (Number.isNaN(value)) ? 0 : value; });
     357                }
     358
     359                md.materials=md.materials.extrude(md);
     360                md.damage=md.damage.extrude(md);
     361                md.mask=md.mask.extrude(md);
     362                md.qmu=md.qmu.extrude(md);
     363                md.basalforcings=md.basalforcings.extrude(md);
     364
     365                //increase connectivity if less than 25:
     366                if (md.mesh.average_vertex_connectivity<=25) {
     367                        md.mesh.average_vertex_connectivity=100;
     368                }
     369        } //}}}
     370        this.extract = function(md,area) { //{{{
     371                //extract - extract a model according to an Argus contour or flag list
     372                //
     373                //   This routine extracts a submodel from a bigger model with respect to a given contour
     374                //   md must be followed by the corresponding exp file or flags list
     375                //   It can either be a domain file (argus type, .exp extension), or an array of element flags.
     376                //   If user wants every element outside the domain to be
     377                //   extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
     378                //   an empty string '' will be considered as an empty domain
     379                //   a string 'all' will be considered as the entire domain
     380                //
     381                //   Usage:
     382                //      md2=extract(md,area);
     383                //
     384                //   Examples:
     385                //      md2=extract(md,'Domain.exp');
     386                //
     387                //   See also: EXTRUDE, COLLAPSE
     388
     389                //some checks on list of arguments
     390                var argc = arguments.length;
     391
     392                if (!((argc == 2) | (argc == 1))) {
     393                        error('extract error message: bad usage');
     394                }
    334395               
    335                                 //md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
    336                                 var indices = ArrayFind(md.mesh.elementconnectivity, 0);
    337                                 for (var i = 0; i < indices.length; i++) {
    338                                         md.mesh.elementconnectivity[i] = NaN;
    339                                 };
     396                //get elements that are inside area
     397                flag_elem=FlagElements(this,area);
     398
     399                if (!ArrayAnyEqual(flag_elem,1)) {
     400                        error('extracted model is empty!');
     401                }
     402
     403                /*kick out all elements with 3 dirichlets: actually, not so fast, not good for javscript usually.
     404                spc_elem=find(~flag_elem);
     405                spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
     406                flag=ones(md1.mesh.numberofvertices,1);
     407                flag(spc_node)=0;
     408                pos=find(sum(flag(md1.mesh.elements),2)==0);
     409                flag_elem(pos)=0;*/
     410
     411                //extracted elements and nodes lists
     412                var pos_elem = ArrayFind(flag_elem,1);
     413                var dup_nodes= new Array(pos_elem.length*3);
     414                for(var i=0;i<pos_elem.length;i++){
     415                        dup_nodes[3*i]=md.mesh.elements[pos_elem[i]][0]-1;
     416                        dup_nodes[3*i+1]=md.mesh.elements[pos_elem[i]][1]-1;
     417                        dup_nodes[3*i+2]=md.mesh.elements[pos_elem[i]][2]-1;
     418                }
     419                pos_node=ArrayUnique(dup_nodes); pos_node=ArraySort(pos_node);
     420
     421                //keep track of some fields
     422                var numberofvertices1=md.mesh.numberofvertices;
     423                var numberofelements1=md.mesh.numberofelements;
     424                var numberofvertices2=pos_node.length;
     425                var numberofelements2=pos_elem.length;
     426                var flag_node=NewArrayFill(numberofvertices1,0);
     427                for (var i=0;i<pos_node.length;i++)flag_node[pos_node[i]]=1;
     428
     429                //Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
     430                Pelem=NewArrayFill(numberofelements1,0);
     431                for (var i=0;i<numberofelements2;i++) Pelem[pos_elem[i]]=i;
     432                Pnode=NewArrayFill(numberofvertices1,0);
     433                for (var i=0;i<numberofvertices2;i++) Pnode[pos_node[i]]=i;
     434                //renumber the elements (some nodes won't exist anymore)
     435                var elements_2=NewArrayFill2D(numberofelements2,3,0);
     436                for (var i=0;i<numberofelements2;i++){
     437                        for (var j=0;j<3;j++){
     438                                elements_2[i][j]=Pnode[md.mesh.elements[i][j]-1]+1;
     439                        }
     440                }
    340441               
    341                                 for (var i = 2; i < numlayers; i++) {
    342                                         //md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
    343                                         //=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
    344                                         for (var j = (i-1)*md.mesh.numberofelements2d; j <= (i)*md.mesh.numberofelements2d-1; j++) {
    345                                                 for (var k = 0; k < 3; k++) {
    346                                                         md.mesh.elementconnectivity[j][k] += md.mesh.numberofelements2d;
     442                /*if isa(md.mesh,'mesh3dprisms'),
     443                        elements_2(:,4)=Pnode(elements_2(:,4));
     444                        elements_2(:,5)=Pnode(elements_2(:,5));
     445                        elements_2(:,6)=Pnode(elements_2(:,6));
     446                end
     447                */
     448               
     449                //OK, now create the new model!
     450
     451                //take every field from model
     452                var md2=md.deepcopy(md);
     453                //var md2=new model(); md2.mesh=new mesh3dsurface();
     454               
     455                //deal with mesh: {{{
     456                md2.mesh.numberofvertices=numberofvertices2;
     457                md2.mesh.numberofelements=numberofelements2;
     458                md2.mesh.elements=elements_2;
     459
     460                md2.mesh.x=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.x[i]=md.mesh.x[pos_node[i]];
     461                md2.mesh.y=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.y[i]=md.mesh.y[pos_node[i]];
     462                md2.mesh.z=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.z[i]=md.mesh.z[pos_node[i]];
     463                md2.mesh.lat=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.lat[i]=md.mesh.lat[pos_node[i]];
     464                md2.mesh.long=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.long[i]=md.mesh.long[pos_node[i]];
     465                md2.mesh.r=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.r[i]=md.mesh.r[pos_node[i]];
     466                //}}}
     467                //deal with geometry: {{{
     468                md2.geometry.base=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.base[i]=md.geometry.base[pos_node[i]];
     469                md2.geometry.thickness=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.thickness[i]=md.geometry.thickness[pos_node[i]];
     470                md2.geometry.surface=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.surface[i]=md.geometry.surface[pos_node[i]];
     471                md2.geometry.bed=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.bed[i]=md.geometry.bed[pos_node[i]];
     472                //}}}
     473
     474                //Keep track of pos_node and pos_elem
     475                for (var i=0;i<md2.mesh.numberofvertices;i++)pos_node[i]=pos_node[i]+1;
     476                for (var i=0;i<md2.mesh.numberofelements;i++)pos_elem[i]=pos_elem[i]+1;
     477                md2.mesh.extractedvertices=pos_node;
     478                md2.mesh.extractedelements=pos_elem;
     479
     480                return md2;
     481
     482                //automatically modify fields
     483
     484                //loop over model fields
     485                //model_fields=fields(md);
     486        } //}}}
     487        this.collapse = function(md) { //{{{
     488                /*
     489                 *COLLAPSE - collapses a 3d mesh into a 2d mesh
     490                 *
     491                 *   This routine collapses a 3d model into a 2d model
     492                 *   and collapses all the fileds of the 3d model by
     493                 *   taking their depth-averaged values
     494                 *
     495                 *   Usage:
     496                 *       md=collapse(md)
     497                 *
     498                 *   See also: EXTRUDE, MODELEXTRACT
     499                 */
     500
     501                // Check that the model is really a 3d model
     502                if (md.mesh.elementtype() !== 'Penta') {
     503                        console.error('collapse error message: only 3d mesh can be collapsed')
     504                }
     505
     506                // Start with changing all the fields from the 3d mesh
     507
     508                // dealing with the friction law
     509                // drag is limited to nodes that are on the bedrock.
     510                if (md.friction.classname() === 'friction') {
     511                        md.friction.coefficient=project2d(md,md.friction.coefficient,1);
     512                        md.friction.p=project2d(md,md.friction.p,1);
     513                        md.friction.q=project2d(md,md.friction.q,1);
     514                } else if (md.friction.classname() === 'frictionhydro') {
     515                        md.friction.q=project2d(md,md.friction.q,1);
     516                        md.friction.C=project2d(md,md.friction.C,1);
     517                        md.friction.As=project2d(md,md.friction.As,1);
     518                        md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
     519                } else if (md.friction.classname() === 'frictionwaterlayer') {
     520                        md.friction.coefficient=project2d(md,md.friction.coefficient,1);
     521                        md.friction.p=project2d(md,md.friction.p,1);
     522                        md.friction.q=project2d(md,md.friction.q,1);
     523                        md.friction.water_layer=project2d(md,md.friction.water_layer,1);
     524                } else if (md.friction.classname() === 'frictionweertman') {
     525                        md.friction.C=project2d(md,md.friction.C,1);
     526                        md.friction.m=project2d(md,md.friction.m,1);
     527                } else if (md.friction.classname() === 'frictionweertmantemp') {
     528                        md.friction.C=project2d(md,md.friction.C,1);
     529                        md.friction.m=project2d(md,md.friction.m,1);
     530                } else {
     531                        disp('friction type not supported');
     532                }
     533
     534                // observations
     535                if (!Number.isNaN(md.inversion.vx_obs))
     536                        md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
     537
     538                if (!Number.isNaN(md.inversion.vy_obs))
     539                        md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
     540
     541                if (!Number.isNaN(md.inversion.vel_obs))
     542                        md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
     543
     544                if (!Number.isNaN(md.inversion.cost_functions_coefficients))
     545                        md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
     546
     547                if (numel(md.inversion.min_parameters)>1)
     548                        md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
     549
     550                if (numel(md.inversion.max_parameters)>1)
     551                        md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
     552
     553                if (md.smb.classname() === 'SMBforcing' && !Number.isNaN(md.smb.mass_balance)) {
     554                        md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
     555                } else if (md.smb.classname() === 'SMBhenning' && !Number.isNaN(md.smb.smbref)) {
     556                        md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
     557                }
     558
     559                // results
     560                if (!Number.isNaN(md.initialization.vx))
     561                        md.initialization.vx=DepthAverage(md,md.initialization.vx);
     562                if (!Number.isNaN(md.initialization.vy))
     563                        md.initialization.vy=DepthAverage(md,md.initialization.vy);
     564                if (!Number.isNaN(md.initialization.vz))
     565                        md.initialization.vz=DepthAverage(md,md.initialization.vz);
     566                if (!Number.isNaN(md.initialization.vel))
     567                        md.initialization.vel=DepthAverage(md,md.initialization.vel);
     568                if (!Number.isNaN(md.initialization.temperature))
     569                        md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
     570                if (!Number.isNaN(md.initialization.pressure))
     571                        md.initialization.pressure=project2d(md,md.initialization.pressure,1);
     572                if (!Number.isNaN(md.initialization.sediment_head))
     573                        md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
     574                if (!Number.isNaN(md.initialization.epl_head))
     575                        md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
     576                if (!Number.isNaN(md.initialization.epl_thickness))
     577                        md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
     578
     579                // giaivins
     580                if (!Number.isNaN(md.gia.mantle_viscosity))
     581                        md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1);
     582                if (!Number.isNaN(md.gia.lithosphere_thickness))
     583                        md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1);
     584
     585                // elementstype
     586                if (!Number.isNaN(md.flowequation.element_equation)) {
     587                        md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
     588                        md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
     589                        md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
     590                        md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
     591                        md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
     592                }
     593
     594                // boundary conditions
     595                md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
     596                md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
     597                md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
     598                md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
     599                md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
     600                md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
     601                if (!Number.isNaN(md.damage.spcdamage))
     602                        md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
     603                md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
     604
     605                // Hydrologydc variables
     606                if (md.hydrology.classname() === 'hydrologydc') {
     607                        md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
     608                        md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
     609                        md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
     610                        md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
     611
     612                        if(md.hydrology.isefficientlayer==1)
     613                                md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
     614                }
     615               
     616                // materials
     617                md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
     618                md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
     619               
     620                // damage:
     621                if (md.damage.isdamage)
     622                        md.damage.D=DepthAverage(md,md.damage.D);
     623
     624                // special for thermal modeling:
     625                if (!Number.isNaN(md.basalforcings.groundedice_melting_rate))
     626                        md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1);
     627
     628                if (!Number.isNaN(md.basalforcings.floatingice_melting_rate))
     629                        md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1);
     630
     631                md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); // bedrock only gets geothermal flux
     632
     633                // update of connectivity matrix
     634                md.mesh.average_vertex_connectivity=25;
     635
     636                // Collapse the mesh
     637                var nodes2d=md.mesh.numberofvertices2d;
     638                var elements2d=md.mesh.numberofelements2d;
     639
     640                // parameters
     641                md.geometry.surface=project2d(md,md.geometry.surface,1);
     642                md.geometry.thickness=project2d(md,md.geometry.thickness,1);
     643                md.geometry.base=project2d(md,md.geometry.base,1);
     644                if (!Number.isNaN(md.geometry.bed))
     645                        md.geometry.bed=project2d(md,md.geometry.bed,1);
     646
     647                if (!Number.isNaN(md.mask.ocean_levelset))
     648                        md.mask.ocean_levelset=project2d(md,md.mask.ocean_levelset,1);
     649
     650                if (!Number.isNaN(md.mask.ice_levelset))
     651                        md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
     652
     653                // lat long
     654                if (numel(md.mesh.lat) === md.mesh.numberofvertices)
     655                        md.mesh.lat=project2d(md,md.mesh.lat,1);
     656                if (numel(md.mesh.long) === md.mesh.numberofvertices)
     657                        md.mesh.long=project2d(md,md.mesh.long,1);
     658                if (numel(md.mesh.scale_factor) === md.mesh.numberofvertices)
     659                        md.mesh.scale_factor=project2d(md,md.mesh.scale_factor,1);
     660
     661                // Initialize with the 2d mesh
     662                var mesh = new mesh2d();
     663                mesh.x=md.mesh.x2d;
     664                mesh.y=md.mesh.y2d;
     665                mesh.numberofvertices=md.mesh.numberofvertices2d;
     666                mesh.numberofelements=md.mesh.numberofelements2d;
     667                mesh.elements=md.mesh.elements2d;
     668
     669                if (!Number.isNaN(md.mesh.vertexonboundary))
     670                        mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
     671                if (!Number.isNaN(md.mesh.elementconnectivity))
     672                        mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
     673
     674                md.mesh=mesh;
     675                md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
     676                md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
     677                md.mesh.segments=contourenvelope(md.mesh);
     678
     679                return md;
     680        } /*}}}*/
     681        this.deepCopy = function(md) { //{{{
     682                /*
     683                 *DEEPCOPY - returns a deep copy of the model.
     684                 *
     685                 *   This routine creates a new model with new objects
     686                 *   for each corresponding object in the original model,
     687                 *   so that changes in one do not affect the other.
     688                 *
     689                 *   Usage:
     690                 *       md1=deepCopy(md)
     691                 *
     692                 */
     693                function recursiveDeepCopy(obj) {
     694                        var returnValue;
     695
     696                        switch (typeof obj) {
     697                                case "object":
     698                                        if (obj === null) {
     699                                                // null => null
     700                                                returnValue = null;
     701                                        } else {
     702                                                switch (toString.call(obj)) {
     703                                                        case "[object Array]":
     704                                                                // It's an array, create a new array with deep copies of the entries
     705                                                                returnValue = obj.map(recursiveDeepCopy);
     706                                                                break;
     707                                                        default:
     708                                                                // Some other kind of object, deep-copy its properties into a new object
     709                                                                returnValue = Object.keys(obj).reduce(function(prev, key) {
     710                                                                        prev[key] = recursiveDeepCopy(obj[key]);
     711                                                                        return prev;
     712                                                                }, {});
     713                                                                break;
    347714                                                }
    348715                                        }
    349                                 }
    350                                 md.mesh.elementconnectivity = md.mesh.elementconnectivity.map(function(value) { return (Number.isNaN(value)) ? 0 : value; });
    351                         }
    352 
    353                         md.materials=md.materials.extrude(md);
    354                         md.damage=md.damage.extrude(md);
    355                         md.mask=md.mask.extrude(md);
    356                         md.qmu=md.qmu.extrude(md);
    357                         md.basalforcings=md.basalforcings.extrude(md);
    358 
    359                         //increase connectivity if less than 25:
    360                         if (md.mesh.average_vertex_connectivity<=25) {
    361                                 md.mesh.average_vertex_connectivity=100;
    362                         }
    363                 } //}}}
    364                 this.extract = function(md,area) { //{{{
    365                         //extract - extract a model according to an Argus contour or flag list
    366                         //
    367                         //   This routine extracts a submodel from a bigger model with respect to a given contour
    368                         //   md must be followed by the corresponding exp file or flags list
    369                         //   It can either be a domain file (argus type, .exp extension), or an array of element flags.
    370                         //   If user wants every element outside the domain to be
    371                         //   extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
    372                         //   an empty string '' will be considered as an empty domain
    373                         //   a string 'all' will be considered as the entire domain
    374                         //
    375                         //   Usage:
    376                         //      md2=extract(md,area);
    377                         //
    378                         //   Examples:
    379                         //      md2=extract(md,'Domain.exp');
    380                         //
    381                         //   See also: EXTRUDE, COLLAPSE
    382 
    383                         //some checks on list of arguments
    384                         var argc = arguments.length;
    385 
    386                         if (!((argc == 2) | (argc == 1))) {
    387                                 throw "extract error message: bad usage";
    388                         }
    389                        
    390                         //get elements that are inside area
    391                         flag_elem=FlagElements(this,area);
    392 
    393                         if (!ArrayAnyEqual(flag_elem,1))throw "extracted model is empty!"
    394 
    395                         /*kick out all elements with 3 dirichlets: actually, not so fast, not good for javscript usually.
    396                         spc_elem=find(~flag_elem);
    397                         spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
    398                         flag=ones(md1.mesh.numberofvertices,1);
    399                         flag(spc_node)=0;
    400                         pos=find(sum(flag(md1.mesh.elements),2)==0);
    401                         flag_elem(pos)=0;*/
    402 
    403                         //extracted elements and nodes lists
    404                         var pos_elem = ArrayFind(flag_elem,1);
    405                         var dup_nodes= new Array(pos_elem.length*3);
    406                         for(var i=0;i<pos_elem.length;i++){
    407                                 dup_nodes[3*i]=md.mesh.elements[pos_elem[i]][0]-1;
    408                                 dup_nodes[3*i+1]=md.mesh.elements[pos_elem[i]][1]-1;
    409                                 dup_nodes[3*i+2]=md.mesh.elements[pos_elem[i]][2]-1;
    410                         }
    411                         pos_node=ArrayUnique(dup_nodes); pos_node=ArraySort(pos_node);
    412 
    413                         //keep track of some fields
    414                         var numberofvertices1=md.mesh.numberofvertices;
    415                         var numberofelements1=md.mesh.numberofelements;
    416                         var numberofvertices2=pos_node.length;
    417                         var numberofelements2=pos_elem.length;
    418                         var flag_node=NewArrayFill(numberofvertices1,0);
    419                         for (var i=0;i<pos_node.length;i++)flag_node[pos_node[i]]=1;
    420 
    421                         //Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
    422                         Pelem=NewArrayFill(numberofelements1,0);
    423                         for (var i=0;i<numberofelements2;i++) Pelem[pos_elem[i]]=i;
    424                         Pnode=NewArrayFill(numberofvertices1,0);
    425                         for (var i=0;i<numberofvertices2;i++) Pnode[pos_node[i]]=i;
    426                         //renumber the elements (some nodes won't exist anymore)
    427                         var elements_2=NewArrayFill2D(numberofelements2,3,0);
    428                         for (var i=0;i<numberofelements2;i++){
    429                                 for (var j=0;j<3;j++){
    430                                         elements_2[i][j]=Pnode[md.mesh.elements[i][j]-1]+1;
    431                                 }
    432                         }
    433                        
    434                         /*if isa(md.mesh,'mesh3dprisms'),
    435                                 elements_2(:,4)=Pnode(elements_2(:,4));
    436                                 elements_2(:,5)=Pnode(elements_2(:,5));
    437                                 elements_2(:,6)=Pnode(elements_2(:,6));
    438                         end
    439                         */
    440                        
    441                         //OK, now create the new model!
    442 
    443                         //take every field from model
    444                         var md2=md.deepcopy(md);
    445                         //var md2=new model(); md2.mesh=new mesh3dsurface();
    446                        
    447                         //deal with mesh: {{{
    448                         md2.mesh.numberofvertices=numberofvertices2;
    449                         md2.mesh.numberofelements=numberofelements2;
    450                         md2.mesh.elements=elements_2;
    451 
    452                         md2.mesh.x=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.x[i]=md.mesh.x[pos_node[i]];
    453                         md2.mesh.y=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.y[i]=md.mesh.y[pos_node[i]];
    454                         md2.mesh.z=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.z[i]=md.mesh.z[pos_node[i]];
    455                         md2.mesh.lat=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.lat[i]=md.mesh.lat[pos_node[i]];
    456                         md2.mesh.long=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.long[i]=md.mesh.long[pos_node[i]];
    457                         md2.mesh.r=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.r[i]=md.mesh.r[pos_node[i]];
    458                         //}}}
    459                         //deal with geometry: {{{
    460                         md2.geometry.base=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.base[i]=md.geometry.base[pos_node[i]];
    461                         md2.geometry.thickness=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.thickness[i]=md.geometry.thickness[pos_node[i]];
    462                         md2.geometry.surface=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.surface[i]=md.geometry.surface[pos_node[i]];
    463                         md2.geometry.bed=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.bed[i]=md.geometry.bed[pos_node[i]];
    464                         //}}}
    465 
    466                         //Keep track of pos_node and pos_elem
    467                         for (var i=0;i<md2.mesh.numberofvertices;i++)pos_node[i]=pos_node[i]+1;
    468                         for (var i=0;i<md2.mesh.numberofelements;i++)pos_elem[i]=pos_elem[i]+1;
    469                         md2.mesh.extractedvertices=pos_node;
    470                         md2.mesh.extractedelements=pos_elem;
    471 
    472                         return md2;
    473 
    474 
    475 
    476                         //automatically modify fields
    477 
    478                         //loop over model fields
    479                         //model_fields=fields(md);
    480 
    481 
    482 
    483 
    484 
    485                 } //}}}
    486                 this.collapse = function(md) { //{{{
    487                         /*
    488                          *COLLAPSE - collapses a 3d mesh into a 2d mesh
    489                          *
    490                          *   This routine collapses a 3d model into a 2d model
    491                          *   and collapses all the fileds of the 3d model by
    492                          *   taking their depth-averaged values
    493                          *
    494                          *   Usage:
    495                          *       md=collapse(md)
    496                          *
    497                          *   See also: EXTRUDE, MODELEXTRACT
    498                          */
    499 
    500                         // Check that the model is really a 3d model
    501                         if (md.mesh.elementtype() !== 'Penta') {
    502                                 console.error('collapse error message: only 3d mesh can be collapsed')
    503                         }
    504 
    505                         // Start with changing all the fields from the 3d mesh
    506 
    507                         // dealing with the friction law
    508                         // drag is limited to nodes that are on the bedrock.
    509                         if (md.friction.classname() === 'friction') {
    510                                 md.friction.coefficient=project2d(md,md.friction.coefficient,1);
    511                                 md.friction.p=project2d(md,md.friction.p,1);
    512                                 md.friction.q=project2d(md,md.friction.q,1);
    513                         } else if (md.friction.classname() === 'frictionhydro') {
    514                                 md.friction.q=project2d(md,md.friction.q,1);
    515                                 md.friction.C=project2d(md,md.friction.C,1);
    516                                 md.friction.As=project2d(md,md.friction.As,1);
    517                                 md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
    518                         } else if (md.friction.classname() === 'frictionwaterlayer') {
    519                                 md.friction.coefficient=project2d(md,md.friction.coefficient,1);
    520                                 md.friction.p=project2d(md,md.friction.p,1);
    521                                 md.friction.q=project2d(md,md.friction.q,1);
    522                                 md.friction.water_layer=project2d(md,md.friction.water_layer,1);
    523                         } else if (md.friction.classname() === 'frictionweertman') {
    524                                 md.friction.C=project2d(md,md.friction.C,1);
    525                                 md.friction.m=project2d(md,md.friction.m,1);
    526                         } else if (md.friction.classname() === 'frictionweertmantemp') {
    527                                 md.friction.C=project2d(md,md.friction.C,1);
    528                                 md.friction.m=project2d(md,md.friction.m,1);
    529                         } else {
    530                                 disp('friction type not supported');
    531                         }
    532 
    533                         // observations
    534                         if (!Number.isNaN(md.inversion.vx_obs))
    535                                 md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
    536 
    537                         if (!Number.isNaN(md.inversion.vy_obs))
    538                                 md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
    539 
    540                         if (!Number.isNaN(md.inversion.vel_obs))
    541                                 md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
    542 
    543                         if (!Number.isNaN(md.inversion.cost_functions_coefficients))
    544                                 md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
    545 
    546                         if (numel(md.inversion.min_parameters)>1)
    547                                 md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
    548 
    549                         if (numel(md.inversion.max_parameters)>1)
    550                                 md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
    551 
    552                         if (md.smb.classname() === 'SMBforcing' && !Number.isNaN(md.smb.mass_balance)) {
    553                                 md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
    554                         } else if (md.smb.classname() === 'SMBhenning' && !Number.isNaN(md.smb.smbref)) {
    555                                 md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
    556                         }
    557 
    558                         // results
    559                         if (!Number.isNaN(md.initialization.vx))
    560                                 md.initialization.vx=DepthAverage(md,md.initialization.vx);
    561                         if (!Number.isNaN(md.initialization.vy))
    562                                 md.initialization.vy=DepthAverage(md,md.initialization.vy);
    563                         if (!Number.isNaN(md.initialization.vz))
    564                                 md.initialization.vz=DepthAverage(md,md.initialization.vz);
    565                         if (!Number.isNaN(md.initialization.vel))
    566                                 md.initialization.vel=DepthAverage(md,md.initialization.vel);
    567                         if (!Number.isNaN(md.initialization.temperature))
    568                                 md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
    569                         if (!Number.isNaN(md.initialization.pressure))
    570                                 md.initialization.pressure=project2d(md,md.initialization.pressure,1);
    571                         if (!Number.isNaN(md.initialization.sediment_head))
    572                                 md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
    573                         if (!Number.isNaN(md.initialization.epl_head))
    574                                 md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
    575                         if (!Number.isNaN(md.initialization.epl_thickness))
    576                                 md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
    577 
    578                         // giaivins
    579                         if (!Number.isNaN(md.gia.mantle_viscosity))
    580                                 md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1);
    581                         if (!Number.isNaN(md.gia.lithosphere_thickness))
    582                                 md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1);
    583 
    584                         // elementstype
    585                         if (!Number.isNaN(md.flowequation.element_equation)) {
    586                                 md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
    587                                 md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
    588                                 md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
    589                                 md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
    590                                 md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
    591                         }
    592 
    593                         // boundary conditions
    594                         md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
    595                         md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
    596                         md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
    597                         md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
    598                         md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
    599                         md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
    600                         if (!Number.isNaN(md.damage.spcdamage))
    601                                 md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
    602                         md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
    603 
    604                         // Hydrologydc variables
    605                         if (md.hydrology.classname() === 'hydrologydc') {
    606                                 md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
    607                                 md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
    608                                 md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
    609                                 md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
    610 
    611                                 if(md.hydrology.isefficientlayer==1)
    612                                         md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
    613                         }
    614                        
    615                         // materials
    616                         md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
    617                         md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
    618                        
    619                         // damage:
    620                         if (md.damage.isdamage)
    621                                 md.damage.D=DepthAverage(md,md.damage.D);
    622 
    623                         // special for thermal modeling:
    624                         if (!Number.isNaN(md.basalforcings.groundedice_melting_rate))
    625                                 md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1);
    626 
    627                         if (!Number.isNaN(md.basalforcings.floatingice_melting_rate))
    628                                 md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_melting_rate,1);
    629 
    630                         md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); // bedrock only gets geothermal flux
    631 
    632                         // update of connectivity matrix
    633                         md.mesh.average_vertex_connectivity=25;
    634 
    635                         // Collapse the mesh
    636                         var nodes2d=md.mesh.numberofvertices2d;
    637                         var elements2d=md.mesh.numberofelements2d;
    638 
    639                         // parameters
    640                         md.geometry.surface=project2d(md,md.geometry.surface,1);
    641                         md.geometry.thickness=project2d(md,md.geometry.thickness,1);
    642                         md.geometry.base=project2d(md,md.geometry.base,1);
    643                         if (!Number.isNaN(md.geometry.bed))
    644                                 md.geometry.bed=project2d(md,md.geometry.bed,1);
    645 
    646                         if (!Number.isNaN(md.mask.ocean_levelset))
    647                                 md.mask.ocean_levelset=project2d(md,md.mask.ocean_levelset,1);
    648 
    649                         if (!Number.isNaN(md.mask.ice_levelset))
    650                                 md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
    651 
    652                         // lat long
    653                         if (numel(md.mesh.lat) === md.mesh.numberofvertices)
    654                                 md.mesh.lat=project2d(md,md.mesh.lat,1);
    655                         if (numel(md.mesh.long) === md.mesh.numberofvertices)
    656                                 md.mesh.long=project2d(md,md.mesh.long,1);
    657                         if (numel(md.mesh.scale_factor) === md.mesh.numberofvertices)
    658                                 md.mesh.scale_factor=project2d(md,md.mesh.scale_factor,1);
    659 
    660                         // Initialize with the 2d mesh
    661                         var mesh = new mesh2d();
    662                         mesh.x=md.mesh.x2d;
    663                         mesh.y=md.mesh.y2d;
    664                         mesh.numberofvertices=md.mesh.numberofvertices2d;
    665                         mesh.numberofelements=md.mesh.numberofelements2d;
    666                         mesh.elements=md.mesh.elements2d;
    667 
    668                         if (!Number.isNaN(md.mesh.vertexonboundary))
    669                                 mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
    670                         if (!Number.isNaN(md.mesh.elementconnectivity))
    671                                 mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
    672 
    673                         md.mesh=mesh;
    674                         md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
    675                         md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
    676                         md.mesh.segments=contourenvelope(md.mesh);
    677 
    678                         return md;
    679                 } /*}}}*/
    680                 this.deepCopy = function(md) { //{{{
    681                         /*
    682                          *DEEPCOPY - returns a deep copy of the model.
    683                          *
    684                          *   This routine creates a new model with new objects
    685                          *   for each corresponding object in the original model,
    686                          *   so that changes in one do not affect the other.
    687                          *
    688                          *   Usage:
    689                          *       md1=deepCopy(md)
    690                          *
    691                          */
    692                         function recursiveDeepCopy(obj) {
    693                                 var returnValue;
    694 
    695                                 switch (typeof obj) {
    696                                         case "object":
    697                                                 if (obj === null) {
    698                                                         // null => null
    699                                                         returnValue = null;
    700                                                 } else {
    701                                                         switch (toString.call(obj)) {
    702                                                                 case "[object Array]":
    703                                                                         // It's an array, create a new array with deep copies of the entries
    704                                                                         returnValue = obj.map(recursiveDeepCopy);
    705                                                                         break;
    706                                                                 default:
    707                                                                         // Some other kind of object, deep-copy its properties into a new object
    708                                                                         returnValue = Object.keys(obj).reduce(function(prev, key) {
    709                                                                                 prev[key] = recursiveDeepCopy(obj[key]);
    710                                                                                 return prev;
    711                                                                         }, {});
    712                                                                         break;
    713                                                         }
    714                                                 }
    715                                                 break;
    716                                         default:
    717                                                 // It's a primitive, copy via assignment
    718                                                 returnValue = obj;
    719                                                 break;
    720                                 }
    721                                 return returnValue;
    722                         }
    723                        
    724                         return recursiveDeepCopy(md);
    725                 } /*}}}*/
    726         //properties
    727         // {{{
    728                 //Careful here: no other class should be used as default value this is a bug of matlab
    729                 this.mesh             = 0;
    730                 this.mask             = 0;
    731 
    732                 this.geometry         = 0;
    733                 this.constants        = 0;
    734                 this.smb              = 0;
    735                 this.basalforcings    = 0;
    736                 this.materials        = 0;
    737                 this.damage           = 0;
    738                 this.friction         = 0;
    739                 this.flowequation     = 0;
    740                 this.timestepping     = 0;
    741                 this.initialization   = 0;
    742                 this.rifts            = 0;
    743                 this.slr              = 0;
    744                 this.dsl              = 0;
    745 
    746                 this.debug            = 0;
    747                 this.verbose          = 0;
    748                 this.settings         = 0;
    749                 this.toolkits         = 0;
    750                 this.cluster          = 0;
    751 
    752                 this.balancethickness = 0;
    753                 this.stressbalance    = 0;
    754                 this.groundingline    = 0;
    755                 this.hydrology        = 0;
    756                 this.masstransport    = 0;
    757                 this.thermal          = 0;
    758                 this.steadystate      = 0;
    759                 this.trans            = 0;
    760                 this.levelset         = 0;
    761                 this.calving          = 0;
    762                 this.frontforcings        = 0;
    763                 this.love             = 0;
    764                 this.gia              = 0;
    765                 this.esa              = 0;
    766 
    767                 this.autodiff         = 0;
    768                 this.inversion        = 0;
    769                 this.qmu              = 0;
    770                 this.amr                                         = 0;
    771 
    772                 this.results          = 0;
    773                 this.outputdefinition = 0;
    774                 this.radaroverlay     = 0;
    775                 this.miscellaneous    = 0;
    776                 this.priv             = 0;
    777 
    778                 //set default values for fields
    779                 this.setdefaultparameters();
    780 
    781                 //}}}
     716                                        break;
     717                                default:
     718                                        // It's a primitive, copy via assignment
     719                                        returnValue = obj;
     720                                        break;
     721                        }
     722                        return returnValue;
     723                }
     724               
     725                return recursiveDeepCopy(md);
     726        } /*}}}*/
     727//properties
     728// {{{
     729        //Careful here: no other class should be used as default value this is a bug of matlab
     730        this.mesh             = 0;
     731        this.mask             = 0;
     732
     733        this.geometry         = 0;
     734        this.constants        = 0;
     735        this.smb              = 0;
     736        this.basalforcings    = 0;
     737        this.materials        = 0;
     738        this.damage           = 0;
     739        this.friction         = 0;
     740        this.flowequation     = 0;
     741        this.timestepping     = 0;
     742        this.initialization   = 0;
     743        this.rifts            = 0;
     744        this.dsl              = 0;
     745        this.solidearth       = 0;
     746
     747        this.debug            = 0;
     748        this.verbose          = 0;
     749        this.settings         = 0;
     750        this.toolkits         = 0;
     751        this.cluster          = 0;
     752
     753        this.balancethickness = 0;
     754        this.stressbalance    = 0;
     755        this.groundingline    = 0;
     756        this.hydrology        = 0;
     757        this.masstransport    = 0;
     758        this.thermal          = 0;
     759        this.steadystate      = 0;
     760        this.transient        = 0;
     761        this.levelset         = 0;
     762        this.calving          = 0;
     763        this.frontalforcings  = 0;
     764        this.love             = 0;
     765        this.esa              = 0;
     766        this.sampling         = 0;
     767
     768        this.autodiff         = 0;
     769        this.inversion        = 0;
     770        this.qmu              = 0;
     771        this.amr              = 0;
     772        this.results          = 0;
     773        this.outputdefinition = 0;
     774        this.radaroverlay     = 0;
     775        this.miscellaneous    = 0;
     776        this.priv             = 0;
     777
     778        // Set default values for fields
     779        if (arguments.length == 0) {
     780                this.setdefaultparameters('earth');
     781        } else {
     782                this.setdefaultparameters(planet);
     783        }
     784        //}}}
    782785}
Note: See TracChangeset for help on using the changeset viewer.