Changeset 18492


Ignore:
Timestamp:
09/10/14 12:40:38 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: added Sea ice model (not working yet) into ISSM

Location:
issm/trunk-jpl
Files:
43 added
24 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/m4/analyses.m4

    r18181 r18492  
    1111
    1212dnl with-AdjointBalancethickness{{{
     13
    1314AC_ARG_WITH([AdjointBalancethickness],
     15
    1416        AS_HELP_STRING([--with-AdjointBalancethickness = YES], [compile with AdjointBalancethickness capabilities (default is yes)]),
     17
    1518        [ADJOINTBALANCETHICKNESS=$withval],[ADJOINTBALANCETHICKNESS=yes])
     19
    1620AC_MSG_CHECKING(for AdjointBalancethickness capability compilation)
    1721
     22
    1823HAVE_ADJOINTBALANCETHICKNESS=no
     24
    1925if test "x$ADJOINTBALANCETHICKNESS" = "xyes"; then
     26
    2027        HAVE_ADJOINTBALANCETHICKNESS=yes
     28
    2129        AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS_],[1],[with AdjointBalancethicknesscapability])
    22 fi
     30
     31fi
     32
    2333AM_CONDITIONAL([ADJOINTBALANCETHICKNESS], [test x$HAVE_ADJOINTBALANCETHICKNESS = xyes])
     34
    2435AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS)
     36
    2537dnl }}}
    2638dnl with-AdjointBalancethickness2{{{
     39
    2740AC_ARG_WITH([AdjointBalancethickness2],
     41
    2842        AS_HELP_STRING([--with-AdjointBalancethickness2 = YES], [compile with AdjointBalancethickness2 capabilities (default is yes)]),
     43
    2944        [ADJOINTBALANCETHICKNESS2=$withval],[ADJOINTBALANCETHICKNESS2=yes])
     45
    3046AC_MSG_CHECKING(for AdjointBalancethickness2 capability compilation)
    3147
     48
    3249HAVE_ADJOINTBALANCETHICKNESS2=no
     50
    3351if test "x$ADJOINTBALANCETHICKNESS2" = "xyes"; then
     52
    3453        HAVE_ADJOINTBALANCETHICKNESS2=yes
     54
    3555        AC_DEFINE([_HAVE_ADJOINTBALANCETHICKNESS2_],[1],[with AdjointBalancethickness2capability])
    36 fi
     56
     57fi
     58
    3759AM_CONDITIONAL([ADJOINTBALANCETHICKNESS2], [test x$HAVE_ADJOINTBALANCETHICKNESS2 = xyes])
     60
    3861AC_MSG_RESULT($HAVE_ADJOINTBALANCETHICKNESS2)
     62
    3963dnl }}}
    4064dnl with-AdjointHoriz{{{
     65
    4166AC_ARG_WITH([AdjointHoriz],
     67
    4268        AS_HELP_STRING([--with-AdjointHoriz = YES], [compile with AdjointHoriz capabilities (default is yes)]),
     69
    4370        [ADJOINTHORIZ=$withval],[ADJOINTHORIZ=yes])
     71
    4472AC_MSG_CHECKING(for AdjointHoriz capability compilation)
    4573
     74
    4675HAVE_ADJOINTHORIZ=no
     76
    4777if test "x$ADJOINTHORIZ" = "xyes"; then
     78
    4879        HAVE_ADJOINTHORIZ=yes
     80
    4981        AC_DEFINE([_HAVE_ADJOINTHORIZ_],[1],[with AdjointHorizcapability])
    50 fi
     82
     83fi
     84
    5185AM_CONDITIONAL([ADJOINTHORIZ], [test x$HAVE_ADJOINTHORIZ = xyes])
     86
    5287AC_MSG_RESULT($HAVE_ADJOINTHORIZ)
     88
    5389dnl }}}
    5490dnl with-Balancethickness{{{
     91
    5592AC_ARG_WITH([Balancethickness],
     93
    5694        AS_HELP_STRING([--with-Balancethickness = YES], [compile with Balancethickness capabilities (default is yes)]),
     95
    5796        [BALANCETHICKNESS=$withval],[BALANCETHICKNESS=yes])
     97
    5898AC_MSG_CHECKING(for Balancethickness capability compilation)
    5999
     100
    60101HAVE_BALANCETHICKNESS=no
     102
    61103if test "x$BALANCETHICKNESS" = "xyes"; then
     104
    62105        HAVE_BALANCETHICKNESS=yes
     106
    63107        AC_DEFINE([_HAVE_BALANCETHICKNESS_],[1],[with Balancethicknesscapability])
    64 fi
     108
     109fi
     110
    65111AM_CONDITIONAL([BALANCETHICKNESS], [test x$HAVE_BALANCETHICKNESS = xyes])
     112
    66113AC_MSG_RESULT($HAVE_BALANCETHICKNESS)
     114
    67115dnl }}}
    68116dnl with-Balancethickness2{{{
     117
    69118AC_ARG_WITH([Balancethickness2],
     119
    70120        AS_HELP_STRING([--with-Balancethickness2 = YES], [compile with Balancethickness2 capabilities (default is yes)]),
     121
    71122        [BALANCETHICKNESS2=$withval],[BALANCETHICKNESS2=yes])
     123
    72124AC_MSG_CHECKING(for Balancethickness2 capability compilation)
    73125
     126
    74127HAVE_BALANCETHICKNESS2=no
     128
    75129if test "x$BALANCETHICKNESS2" = "xyes"; then
     130
    76131        HAVE_BALANCETHICKNESS2=yes
     132
    77133        AC_DEFINE([_HAVE_BALANCETHICKNESS2_],[1],[with Balancethickness2capability])
    78 fi
     134
     135fi
     136
    79137AM_CONDITIONAL([BALANCETHICKNESS2], [test x$HAVE_BALANCETHICKNESS2 = xyes])
     138
    80139AC_MSG_RESULT($HAVE_BALANCETHICKNESS2)
     140
    81141dnl }}}
    82142dnl with-BalancethicknessSoft{{{
     143
    83144AC_ARG_WITH([BalancethicknessSoft],
     145
    84146        AS_HELP_STRING([--with-BalancethicknessSoft = YES], [compile with BalancethicknessSoft capabilities (default is yes)]),
     147
    85148        [BALANCETHICKNESSSOFT=$withval],[BALANCETHICKNESSSOFT=yes])
     149
    86150AC_MSG_CHECKING(for BalancethicknessSoft capability compilation)
    87151
     152
    88153HAVE_BALANCETHICKNESSSOFT=no
     154
    89155if test "x$BALANCETHICKNESSSOFT" = "xyes"; then
     156
    90157        HAVE_BALANCETHICKNESSSOFT=yes
     158
    91159        AC_DEFINE([_HAVE_BALANCETHICKNESSSOFT_],[1],[with BalancethicknessSoftcapability])
    92 fi
     160
     161fi
     162
    93163AM_CONDITIONAL([BALANCETHICKNESSSOFT], [test x$HAVE_BALANCETHICKNESSSOFT = xyes])
     164
    94165AC_MSG_RESULT($HAVE_BALANCETHICKNESSSOFT)
     166
    95167dnl }}}
    96168dnl with-Balancevelocity{{{
     169
    97170AC_ARG_WITH([Balancevelocity],
     171
    98172        AS_HELP_STRING([--with-Balancevelocity = YES], [compile with Balancevelocity capabilities (default is yes)]),
     173
    99174        [BALANCEVELOCITY=$withval],[BALANCEVELOCITY=yes])
     175
    100176AC_MSG_CHECKING(for Balancevelocity capability compilation)
    101177
     178
    102179HAVE_BALANCEVELOCITY=no
     180
    103181if test "x$BALANCEVELOCITY" = "xyes"; then
     182
    104183        HAVE_BALANCEVELOCITY=yes
     184
    105185        AC_DEFINE([_HAVE_BALANCEVELOCITY_],[1],[with Balancevelocitycapability])
    106 fi
     186
     187fi
     188
    107189AM_CONDITIONAL([BALANCEVELOCITY], [test x$HAVE_BALANCEVELOCITY = xyes])
     190
    108191AC_MSG_RESULT($HAVE_BALANCEVELOCITY)
     192
    109193dnl }}}
    110194dnl with-L2ProjectionEPL{{{
     195
    111196AC_ARG_WITH([L2ProjectionEPL],
     197
    112198        AS_HELP_STRING([--with-L2ProjectionEPL = YES], [compile with L2ProjectionEPL capabilities (default is yes)]),
     199
    113200        [L2PROJECTIONEPL=$withval],[L2PROJECTIONEPL=yes])
     201
    114202AC_MSG_CHECKING(for L2ProjectionEPL capability compilation)
    115203
     204
    116205HAVE_L2PROJECTIONEPL=no
     206
    117207if test "x$L2PROJECTIONEPL" = "xyes"; then
     208
    118209        HAVE_L2PROJECTIONEPL=yes
     210
    119211        AC_DEFINE([_HAVE_L2PROJECTIONEPL_],[1],[with L2ProjectionEPLcapability])
    120 fi
     212
     213fi
     214
    121215AM_CONDITIONAL([L2PROJECTIONEPL], [test x$HAVE_L2PROJECTIONEPL = xyes])
     216
    122217AC_MSG_RESULT($HAVE_L2PROJECTIONEPL)
     218
    123219dnl }}}
    124220dnl with-L2ProjectionBase{{{
     221
    125222AC_ARG_WITH([L2ProjectionBase],
     223
    126224        AS_HELP_STRING([--with-L2ProjectionBase = YES], [compile with L2ProjectionBase capabilities (default is yes)]),
     225
    127226        [L2PROJECTIONBASE=$withval],[L2PROJECTIONBASE=yes])
     227
    128228AC_MSG_CHECKING(for L2ProjectionBase capability compilation)
    129229
     230
    130231HAVE_L2PROJECTIONBASE=no
     232
    131233if test "x$L2PROJECTIONBASE" = "xyes"; then
     234
    132235        HAVE_L2PROJECTIONBASE=yes
     236
    133237        AC_DEFINE([_HAVE_L2PROJECTIONBASE_],[1],[with L2ProjectionBasecapability])
    134 fi
     238
     239fi
     240
    135241AM_CONDITIONAL([L2PROJECTIONBASE], [test x$HAVE_L2PROJECTIONBASE = xyes])
     242
    136243AC_MSG_RESULT($HAVE_L2PROJECTIONBASE)
     244
    137245dnl }}}
    138246dnl with-DamageEvolution{{{
     247
    139248AC_ARG_WITH([DamageEvolution],
     249
    140250        AS_HELP_STRING([--with-DamageEvolution = YES], [compile with DamageEvolution capabilities (default is yes)]),
     251
    141252        [DAMAGEEVOLUTION=$withval],[DAMAGEEVOLUTION=yes])
     253
    142254AC_MSG_CHECKING(for DamageEvolution capability compilation)
    143255
     256
    144257HAVE_DAMAGEEVOLUTION=no
     258
    145259if test "x$DAMAGEEVOLUTION" = "xyes"; then
     260
    146261        HAVE_DAMAGEEVOLUTION=yes
     262
    147263        AC_DEFINE([_HAVE_DAMAGEEVOLUTION_],[1],[with DamageEvolutioncapability])
    148 fi
     264
     265fi
     266
    149267AM_CONDITIONAL([DAMAGEEVOLUTION], [test x$HAVE_DAMAGEEVOLUTION = xyes])
     268
    150269AC_MSG_RESULT($HAVE_DAMAGEEVOLUTION)
     270
    151271dnl }}}
    152272dnl with-Stressbalance{{{
     273
    153274AC_ARG_WITH([Stressbalance],
     275
    154276        AS_HELP_STRING([--with-Stressbalance = YES], [compile with Stressbalance capabilities (default is yes)]),
     277
    155278        [STRESSBALANCE=$withval],[STRESSBALANCE=yes])
     279
    156280AC_MSG_CHECKING(for Stressbalance capability compilation)
    157281
     282
    158283HAVE_STRESSBALANCE=no
     284
    159285if test "x$STRESSBALANCE" = "xyes"; then
     286
    160287        HAVE_STRESSBALANCE=yes
     288
    161289        AC_DEFINE([_HAVE_STRESSBALANCE_],[1],[with Stressbalancecapability])
    162 fi
     290
     291fi
     292
    163293AM_CONDITIONAL([STRESSBALANCE], [test x$HAVE_STRESSBALANCE = xyes])
     294
    164295AC_MSG_RESULT($HAVE_STRESSBALANCE)
     296
    165297dnl }}}
    166298dnl with-StressbalanceSIA{{{
     299
    167300AC_ARG_WITH([StressbalanceSIA],
     301
    168302        AS_HELP_STRING([--with-StressbalanceSIA = YES], [compile with StressbalanceSIA capabilities (default is yes)]),
     303
    169304        [STRESSBALANCESIA=$withval],[STRESSBALANCESIA=yes])
     305
    170306AC_MSG_CHECKING(for StressbalanceSIA capability compilation)
    171307
     308
    172309HAVE_STRESSBALANCESIA=no
     310
    173311if test "x$STRESSBALANCESIA" = "xyes"; then
     312
    174313        HAVE_STRESSBALANCESIA=yes
     314
    175315        AC_DEFINE([_HAVE_STRESSBALANCESIA_],[1],[with StressbalanceSIAcapability])
    176 fi
     316
     317fi
     318
    177319AM_CONDITIONAL([STRESSBALANCESIA], [test x$HAVE_STRESSBALANCESIA = xyes])
     320
    178321AC_MSG_RESULT($HAVE_STRESSBALANCESIA)
     322
    179323dnl }}}
    180324dnl with-StressbalanceVertical{{{
     325
    181326AC_ARG_WITH([StressbalanceVertical],
     327
    182328        AS_HELP_STRING([--with-StressbalanceVertical = YES], [compile with StressbalanceVertical capabilities (default is yes)]),
     329
    183330        [STRESSBALANCEVERTICAL=$withval],[STRESSBALANCEVERTICAL=yes])
     331
    184332AC_MSG_CHECKING(for StressbalanceVertical capability compilation)
    185333
     334
    186335HAVE_STRESSBALANCEVERTICAL=no
     336
    187337if test "x$STRESSBALANCEVERTICAL" = "xyes"; then
     338
    188339        HAVE_STRESSBALANCEVERTICAL=yes
     340
    189341        AC_DEFINE([_HAVE_STRESSBALANCEVERTICAL_],[1],[with StressbalanceVerticalcapability])
    190 fi
     342
     343fi
     344
    191345AM_CONDITIONAL([STRESSBALANCEVERTICAL], [test x$HAVE_STRESSBALANCEVERTICAL = xyes])
     346
    192347AC_MSG_RESULT($HAVE_STRESSBALANCEVERTICAL)
     348
    193349dnl }}}
    194350dnl with-Enthalpy{{{
     351
    195352AC_ARG_WITH([Enthalpy],
     353
    196354        AS_HELP_STRING([--with-Enthalpy = YES], [compile with Enthalpy capabilities (default is yes)]),
     355
    197356        [ENTHALPY=$withval],[ENTHALPY=yes])
     357
    198358AC_MSG_CHECKING(for Enthalpy capability compilation)
    199359
     360
    200361HAVE_ENTHALPY=no
     362
    201363if test "x$ENTHALPY" = "xyes"; then
     364
    202365        HAVE_ENTHALPY=yes
     366
    203367        AC_DEFINE([_HAVE_ENTHALPY_],[1],[with Enthalpycapability])
    204 fi
     368
     369fi
     370
    205371AM_CONDITIONAL([ENTHALPY], [test x$HAVE_ENTHALPY = xyes])
     372
    206373AC_MSG_RESULT($HAVE_ENTHALPY)
     374
    207375dnl }}}
    208376dnl with-HydrologyShreve{{{
     377
    209378AC_ARG_WITH([HydrologyShreve],
     379
    210380        AS_HELP_STRING([--with-HydrologyShreve = YES], [compile with HydrologyShreve capabilities (default is yes)]),
     381
    211382        [HYDROLOGYSHREVE=$withval],[HYDROLOGYSHREVE=yes])
     383
    212384AC_MSG_CHECKING(for HydrologyShreve capability compilation)
    213385
     386
    214387HAVE_HYDROLOGYSHREVE=no
     388
    215389if test "x$HYDROLOGYSHREVE" = "xyes"; then
     390
    216391        HAVE_HYDROLOGYSHREVE=yes
     392
    217393        AC_DEFINE([_HAVE_HYDROLOGYSHREVE_],[1],[with HydrologyShrevecapability])
    218 fi
     394
     395fi
     396
    219397AM_CONDITIONAL([HYDROLOGYSHREVE], [test x$HAVE_HYDROLOGYSHREVE = xyes])
     398
    220399AC_MSG_RESULT($HAVE_HYDROLOGYSHREVE)
     400
    221401dnl }}}
    222402dnl with-HydrologyDCInefficient{{{
     403
    223404AC_ARG_WITH([HydrologyDCInefficient],
     405
    224406        AS_HELP_STRING([--with-HydrologyDCInefficient = YES], [compile with HydrologyDCInefficient capabilities (default is yes)]),
     407
    225408        [HYDROLOGYDCINEFFICIENT=$withval],[HYDROLOGYDCINEFFICIENT=yes])
     409
    226410AC_MSG_CHECKING(for HydrologyDCInefficient capability compilation)
    227411
     412
    228413HAVE_HYDROLOGYDCINEFFICIENT=no
     414
    229415if test "x$HYDROLOGYDCINEFFICIENT" = "xyes"; then
     416
    230417        HAVE_HYDROLOGYDCINEFFICIENT=yes
     418
    231419        AC_DEFINE([_HAVE_HYDROLOGYDCINEFFICIENT_],[1],[with HydrologyDCInefficientcapability])
    232 fi
     420
     421fi
     422
    233423AM_CONDITIONAL([HYDROLOGYDCINEFFICIENT], [test x$HAVE_HYDROLOGYDCINEFFICIENT = xyes])
     424
    234425AC_MSG_RESULT($HAVE_HYDROLOGYDCINEFFICIENT)
     426
    235427dnl }}}
    236428dnl with-HydrologyDCEfficient{{{
     429
    237430AC_ARG_WITH([HydrologyDCEfficient],
     431
    238432        AS_HELP_STRING([--with-HydrologyDCEfficient = YES], [compile with HydrologyDCEfficient capabilities (default is yes)]),
     433
    239434        [HYDROLOGYDCEFFICIENT=$withval],[HYDROLOGYDCEFFICIENT=yes])
     435
    240436AC_MSG_CHECKING(for HydrologyDCEfficient capability compilation)
    241437
     438
    242439HAVE_HYDROLOGYDCEFFICIENT=no
     440
    243441if test "x$HYDROLOGYDCEFFICIENT" = "xyes"; then
     442
    244443        HAVE_HYDROLOGYDCEFFICIENT=yes
     444
    245445        AC_DEFINE([_HAVE_HYDROLOGYDCEFFICIENT_],[1],[with HydrologyDCEfficientcapability])
    246 fi
     446
     447fi
     448
    247449AM_CONDITIONAL([HYDROLOGYDCEFFICIENT], [test x$HAVE_HYDROLOGYDCEFFICIENT = xyes])
     450
    248451AC_MSG_RESULT($HAVE_HYDROLOGYDCEFFICIENT)
     452
    249453dnl }}}
    250454dnl with-Melting{{{
     455
    251456AC_ARG_WITH([Melting],
     457
    252458        AS_HELP_STRING([--with-Melting = YES], [compile with Melting capabilities (default is yes)]),
     459
    253460        [MELTING=$withval],[MELTING=yes])
     461
    254462AC_MSG_CHECKING(for Melting capability compilation)
    255463
     464
    256465HAVE_MELTING=no
     466
    257467if test "x$MELTING" = "xyes"; then
     468
    258469        HAVE_MELTING=yes
     470
    259471        AC_DEFINE([_HAVE_MELTING_],[1],[with Meltingcapability])
    260 fi
     472
     473fi
     474
    261475AM_CONDITIONAL([MELTING], [test x$HAVE_MELTING = xyes])
     476
    262477AC_MSG_RESULT($HAVE_MELTING)
     478
    263479dnl }}}
    264480dnl with-Masstransport{{{
     481
    265482AC_ARG_WITH([Masstransport],
     483
    266484        AS_HELP_STRING([--with-Masstransport = YES], [compile with Masstransport capabilities (default is yes)]),
     485
    267486        [MASSTRANSPORT=$withval],[MASSTRANSPORT=yes])
     487
    268488AC_MSG_CHECKING(for Masstransport capability compilation)
    269489
     490
    270491HAVE_MASSTRANSPORT=no
     492
    271493if test "x$MASSTRANSPORT" = "xyes"; then
     494
    272495        HAVE_MASSTRANSPORT=yes
     496
    273497        AC_DEFINE([_HAVE_MASSTRANSPORT_],[1],[with Masstransportcapability])
    274 fi
     498
     499fi
     500
    275501AM_CONDITIONAL([MASSTRANSPORT], [test x$HAVE_MASSTRANSPORT = xyes])
     502
    276503AC_MSG_RESULT($HAVE_MASSTRANSPORT)
     504
    277505dnl }}}
    278506dnl with-FreeSurfaceBase{{{
     507
    279508AC_ARG_WITH([FreeSurfaceBase],
     509
    280510        AS_HELP_STRING([--with-FreeSurfaceBase = YES], [compile with FreeSurfaceBase capabilities (default is yes)]),
     511
    281512        [FREESURFACEBASE=$withval],[FREESURFACEBASE=yes])
     513
    282514AC_MSG_CHECKING(for FreeSurfaceBase capability compilation)
    283515
     516
    284517HAVE_FREESURFACEBASE=no
     518
    285519if test "x$FREESURFACEBASE" = "xyes"; then
     520
    286521        HAVE_FREESURFACEBASE=yes
     522
    287523        AC_DEFINE([_HAVE_FREESURFACEBASE_],[1],[with FreeSurfaceBasecapability])
    288 fi
     524
     525fi
     526
    289527AM_CONDITIONAL([FREESURFACEBASE], [test x$HAVE_FREESURFACEBASE = xyes])
     528
    290529AC_MSG_RESULT($HAVE_FREESURFACEBASE)
     530
    291531dnl }}}
    292532dnl with-FreeSurfaceTop{{{
     533
    293534AC_ARG_WITH([FreeSurfaceTop],
     535
    294536        AS_HELP_STRING([--with-FreeSurfaceTop = YES], [compile with FreeSurfaceTop capabilities (default is yes)]),
     537
    295538        [FREESURFACETOP=$withval],[FREESURFACETOP=yes])
     539
    296540AC_MSG_CHECKING(for FreeSurfaceTop capability compilation)
    297541
     542
    298543HAVE_FREESURFACETOP=no
     544
    299545if test "x$FREESURFACETOP" = "xyes"; then
     546
    300547        HAVE_FREESURFACETOP=yes
     548
    301549        AC_DEFINE([_HAVE_FREESURFACETOP_],[1],[with FreeSurfaceTopcapability])
    302 fi
     550
     551fi
     552
    303553AM_CONDITIONAL([FREESURFACETOP], [test x$HAVE_FREESURFACETOP = xyes])
     554
    304555AC_MSG_RESULT($HAVE_FREESURFACETOP)
     556
    305557dnl }}}
    306558dnl with-ExtrudeFromBase{{{
     559
    307560AC_ARG_WITH([ExtrudeFromBase],
     561
    308562        AS_HELP_STRING([--with-ExtrudeFromBase = YES], [compile with ExtrudeFromBase capabilities (default is yes)]),
     563
    309564        [EXTRUDEFROMBASE=$withval],[EXTRUDEFROMBASE=yes])
     565
    310566AC_MSG_CHECKING(for ExtrudeFromBase capability compilation)
    311567
     568
    312569HAVE_EXTRUDEFROMBASE=no
     570
    313571if test "x$EXTRUDEFROMBASE" = "xyes"; then
     572
    314573        HAVE_EXTRUDEFROMBASE=yes
     574
    315575        AC_DEFINE([_HAVE_EXTRUDEFROMBASE_],[1],[with ExtrudeFromBasecapability])
    316 fi
     576
     577fi
     578
    317579AM_CONDITIONAL([EXTRUDEFROMBASE], [test x$HAVE_EXTRUDEFROMBASE = xyes])
     580
    318581AC_MSG_RESULT($HAVE_EXTRUDEFROMBASE)
     582
    319583dnl }}}
    320584dnl with-ExtrudeFromTop{{{
     585
    321586AC_ARG_WITH([ExtrudeFromTop],
     587
    322588        AS_HELP_STRING([--with-ExtrudeFromTop = YES], [compile with ExtrudeFromTop capabilities (default is yes)]),
     589
    323590        [EXTRUDEFROMTOP=$withval],[EXTRUDEFROMTOP=yes])
     591
    324592AC_MSG_CHECKING(for ExtrudeFromTop capability compilation)
    325593
     594
    326595HAVE_EXTRUDEFROMTOP=no
     596
    327597if test "x$EXTRUDEFROMTOP" = "xyes"; then
     598
    328599        HAVE_EXTRUDEFROMTOP=yes
     600
    329601        AC_DEFINE([_HAVE_EXTRUDEFROMTOP_],[1],[with ExtrudeFromTopcapability])
    330 fi
     602
     603fi
     604
    331605AM_CONDITIONAL([EXTRUDEFROMTOP], [test x$HAVE_EXTRUDEFROMTOP = xyes])
     606
    332607AC_MSG_RESULT($HAVE_EXTRUDEFROMTOP)
     608
    333609dnl }}}
    334610dnl with-DepthAverage{{{
     611
    335612AC_ARG_WITH([DepthAverage],
     613
    336614        AS_HELP_STRING([--with-DepthAverage = YES], [compile with DepthAverage capabilities (default is yes)]),
     615
    337616        [DEPTHAVERAGE=$withval],[DEPTHAVERAGE=yes])
     617
    338618AC_MSG_CHECKING(for DepthAverage capability compilation)
    339619
     620
    340621HAVE_DEPTHAVERAGE=no
     622
    341623if test "x$DEPTHAVERAGE" = "xyes"; then
     624
    342625        HAVE_DEPTHAVERAGE=yes
     626
    343627        AC_DEFINE([_HAVE_DEPTHAVERAGE_],[1],[with DepthAveragecapability])
    344 fi
     628
     629fi
     630
    345631AM_CONDITIONAL([DEPTHAVERAGE], [test x$HAVE_DEPTHAVERAGE = xyes])
     632
    346633AC_MSG_RESULT($HAVE_DEPTHAVERAGE)
     634
    347635dnl }}}
    348636dnl with-SmoothedSurfaceSlopeX{{{
     637
    349638AC_ARG_WITH([SmoothedSurfaceSlopeX],
     639
    350640        AS_HELP_STRING([--with-SmoothedSurfaceSlopeX = YES], [compile with SmoothedSurfaceSlopeX capabilities (default is yes)]),
     641
    351642        [SMOOTHEDSURFACESLOPEX=$withval],[SMOOTHEDSURFACESLOPEX=yes])
     643
    352644AC_MSG_CHECKING(for SmoothedSurfaceSlopeX capability compilation)
    353645
     646
    354647HAVE_SMOOTHEDSURFACESLOPEX=no
     648
    355649if test "x$SMOOTHEDSURFACESLOPEX" = "xyes"; then
     650
    356651        HAVE_SMOOTHEDSURFACESLOPEX=yes
     652
    357653        AC_DEFINE([_HAVE_SMOOTHEDSURFACESLOPEX_],[1],[with SmoothedSurfaceSlopeXcapability])
    358 fi
     654
     655fi
     656
    359657AM_CONDITIONAL([SMOOTHEDSURFACESLOPEX], [test x$HAVE_SMOOTHEDSURFACESLOPEX = xyes])
     658
    360659AC_MSG_RESULT($HAVE_SMOOTHEDSURFACESLOPEX)
     660
    361661dnl }}}
    362662dnl with-SmoothedSurfaceSlopeY{{{
     663
    363664AC_ARG_WITH([SmoothedSurfaceSlopeY],
     665
    364666        AS_HELP_STRING([--with-SmoothedSurfaceSlopeY = YES], [compile with SmoothedSurfaceSlopeY capabilities (default is yes)]),
     667
    365668        [SMOOTHEDSURFACESLOPEY=$withval],[SMOOTHEDSURFACESLOPEY=yes])
     669
    366670AC_MSG_CHECKING(for SmoothedSurfaceSlopeY capability compilation)
    367671
     672
    368673HAVE_SMOOTHEDSURFACESLOPEY=no
     674
    369675if test "x$SMOOTHEDSURFACESLOPEY" = "xyes"; then
     676
    370677        HAVE_SMOOTHEDSURFACESLOPEY=yes
     678
    371679        AC_DEFINE([_HAVE_SMOOTHEDSURFACESLOPEY_],[1],[with SmoothedSurfaceSlopeYcapability])
    372 fi
     680
     681fi
     682
    373683AM_CONDITIONAL([SMOOTHEDSURFACESLOPEY], [test x$HAVE_SMOOTHEDSURFACESLOPEY = xyes])
     684
    374685AC_MSG_RESULT($HAVE_SMOOTHEDSURFACESLOPEY)
     686
     687dnl }}}
     688dnl with-Thermal{{{
     689
     690AC_ARG_WITH([Thermal],
     691
     692        AS_HELP_STRING([--with-Thermal = YES], [compile with Thermal capabilities (default is yes)]),
     693
     694        [THERMAL=$withval],[THERMAL=yes])
     695
     696AC_MSG_CHECKING(for Thermal capability compilation)
     697
     698
     699HAVE_THERMAL=no
     700
     701if test "x$THERMAL" = "xyes"; then
     702
     703        HAVE_THERMAL=yes
     704
     705        AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermalcapability])
     706
     707fi
     708
     709AM_CONDITIONAL([THERMAL], [test x$HAVE_THERMAL = xyes])
     710
     711AC_MSG_RESULT($HAVE_THERMAL)
     712
    375713dnl }}}
    376714dnl with-UzawaPressure{{{
     715
    377716AC_ARG_WITH([UzawaPressure],
     717
    378718        AS_HELP_STRING([--with-UzawaPressure = YES], [compile with UzawaPressure capabilities (default is yes)]),
     719
    379720        [UZAWAPRESSURE=$withval],[UZAWAPRESSURE=yes])
     721
    380722AC_MSG_CHECKING(for UzawaPressure capability compilation)
    381723
     724
    382725HAVE_UZAWAPRESSURE=no
     726
    383727if test "x$UZAWAPRESSURE" = "xyes"; then
     728
    384729        HAVE_UZAWAPRESSURE=yes
     730
    385731        AC_DEFINE([_HAVE_UZAWAPRESSURE_],[1],[with UzawaPressurecapability])
    386 fi
     732
     733fi
     734
    387735AM_CONDITIONAL([UZAWAPRESSURE], [test x$HAVE_UZAWAPRESSURE = xyes])
     736
    388737AC_MSG_RESULT($HAVE_UZAWAPRESSURE)
    389 dnl }}}
    390 dnl with-Thermal{{{
    391 AC_ARG_WITH([Thermal],
    392         AS_HELP_STRING([--with-Thermal = YES], [compile with Thermal capabilities (default is yes)]),
    393         [THERMAL=$withval],[THERMAL=yes])
    394 AC_MSG_CHECKING(for Thermal capability compilation)
    395 
    396 HAVE_THERMAL=no
    397 if test "x$THERMAL" = "xyes"; then
    398         HAVE_THERMAL=yes
    399         AC_DEFINE([_HAVE_THERMAL_],[1],[with Thermalcapability])
    400 fi
    401 AM_CONDITIONAL([THERMAL], [test x$HAVE_THERMAL = xyes])
    402 AC_MSG_RESULT($HAVE_THERMAL)
     738
    403739dnl }}}
    404740dnl with-Gia{{{
     741
    405742AC_ARG_WITH([Gia],
     743
    406744        AS_HELP_STRING([--with-Gia = YES], [compile with Gia capabilities (default is yes)]),
     745
    407746        [GIA=$withval],[GIA=yes])
     747
    408748AC_MSG_CHECKING(for Gia capability compilation)
    409749
     750
    410751HAVE_GIA=no
     752
    411753if test "x$GIA" = "xyes"; then
     754
    412755        HAVE_GIA=yes
     756
    413757        AC_DEFINE([_HAVE_GIA_],[1],[with Giacapability])
    414 fi
     758
     759fi
     760
    415761AM_CONDITIONAL([GIA], [test x$HAVE_GIA = xyes])
     762
    416763AC_MSG_RESULT($HAVE_GIA)
     764
     765dnl }}}
     766dnl with-Seaice{{{
     767
     768AC_ARG_WITH([Seaice],
     769
     770        AS_HELP_STRING([--with-Seaice = YES], [compile with Seaice capabilities (default is yes)]),
     771
     772        [SEAICE=$withval],[SEAICE=yes])
     773
     774AC_MSG_CHECKING(for Seaice capability compilation)
     775
     776
     777HAVE_SEAICE=no
     778
     779if test "x$SEAICE" = "xyes"; then
     780
     781        HAVE_SEAICE=yes
     782
     783        AC_DEFINE([_HAVE_SEAICE_],[1],[with Seaicecapability])
     784
     785fi
     786
     787AM_CONDITIONAL([SEAICE], [test x$HAVE_SEAICE = xyes])
     788
     789AC_MSG_RESULT($HAVE_SEAICE)
     790
    417791dnl }}}
    418792dnl with-Meshdeformation{{{
     793
    419794AC_ARG_WITH([Meshdeformation],
     795
    420796        AS_HELP_STRING([--with-Meshdeformation = YES], [compile with Meshdeformation capabilities (default is yes)]),
     797
    421798        [MESHDEFORMATION=$withval],[MESHDEFORMATION=yes])
     799
    422800AC_MSG_CHECKING(for Meshdeformation capability compilation)
    423801
     802
    424803HAVE_MESHDEFORMATION=no
     804
    425805if test "x$MESHDEFORMATION" = "xyes"; then
     806
    426807        HAVE_MESHDEFORMATION=yes
     808
    427809        AC_DEFINE([_HAVE_MESHDEFORMATION_],[1],[with Meshdeformationcapability])
    428 fi
     810
     811fi
     812
    429813AM_CONDITIONAL([MESHDEFORMATION], [test x$HAVE_MESHDEFORMATION = xyes])
     814
    430815AC_MSG_RESULT($HAVE_MESHDEFORMATION)
     816
    431817dnl }}}
    432818dnl with-Levelset{{{
     819
    433820AC_ARG_WITH([Levelset],
     821
    434822        AS_HELP_STRING([--with-Levelset = YES], [compile with Levelset capabilities (default is yes)]),
     823
    435824        [LEVELSET=$withval],[LEVELSET=yes])
     825
    436826AC_MSG_CHECKING(for Levelset capability compilation)
    437827
     828
    438829HAVE_LEVELSET=no
     830
    439831if test "x$LEVELSET" = "xyes"; then
     832
    440833        HAVE_LEVELSET=yes
     834
    441835        AC_DEFINE([_HAVE_LEVELSET_],[1],[with Levelsetcapability])
    442 fi
     836
     837fi
     838
    443839AM_CONDITIONAL([LEVELSET], [test x$HAVE_LEVELSET = xyes])
     840
    444841AC_MSG_RESULT($HAVE_LEVELSET)
     842
    445843dnl }}}
    446844dnl with-Extrapolation{{{
     845
    447846AC_ARG_WITH([Extrapolation],
     847
    448848        AS_HELP_STRING([--with-Extrapolation = YES], [compile with Extrapolation capabilities (default is yes)]),
     849
    449850        [EXTRAPOLATION=$withval],[EXTRAPOLATION=yes])
     851
    450852AC_MSG_CHECKING(for Extrapolation capability compilation)
    451853
     854
    452855HAVE_EXTRAPOLATION=no
     856
    453857if test "x$EXTRAPOLATION" = "xyes"; then
     858
    454859        HAVE_EXTRAPOLATION=yes
     860
    455861        AC_DEFINE([_HAVE_EXTRAPOLATION_],[1],[with Extrapolationcapability])
    456 fi
     862
     863fi
     864
    457865AM_CONDITIONAL([EXTRAPOLATION], [test x$HAVE_EXTRAPOLATION = xyes])
     866
    458867AC_MSG_RESULT($HAVE_EXTRAPOLATION)
     868
    459869dnl }}}
    460870dnl with-LsfReinitialization{{{
     871
    461872AC_ARG_WITH([LsfReinitialization],
     873
    462874        AS_HELP_STRING([--with-LsfReinitialization = YES], [compile with LsfReinitialization capabilities (default is yes)]),
     875
    463876        [LSFREINITIALIZATION=$withval],[LSFREINITIALIZATION=yes])
     877
    464878AC_MSG_CHECKING(for LsfReinitialization capability compilation)
    465879
     880
    466881HAVE_LSFREINITIALIZATION=no
     882
    467883if test "x$LSFREINITIALIZATION" = "xyes"; then
     884
    468885        HAVE_LSFREINITIALIZATION=yes
     886
    469887        AC_DEFINE([_HAVE_LSFREINITIALIZATION_],[1],[with LsfReinitializationcapability])
    470 fi
     888
     889fi
     890
    471891AM_CONDITIONAL([LSFREINITIALIZATION], [test x$HAVE_LSFREINITIALIZATION = xyes])
     892
    472893AC_MSG_RESULT($HAVE_LSFREINITIALIZATION)
     894
    473895dnl }}}
    474896

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

    r18370 r18492  
    108108                                        ./classes/Materials/Matice.h\
    109109                                        ./classes/Materials/Matice.cpp\
     110                                        ./classes/Materials/Matseaice.h\
     111                                        ./classes/Materials/Matseaice.cpp\
    110112                                        ./classes/Materials/Matpar.h\
    111113                                        ./classes/Materials/Matpar.cpp\
    … …  
    649651issm_sources += ./analyses/DepthAverageAnalysis.cpp
    650652endif
     653if SEAICE
     654issm_sources += ./cores/seaice_core.cpp
     655issm_sources += ./analyses/SeaiceAnalysis.cpp
     656endif
    651657if THERMAL
    652658issm_sources += ./analyses/ThermalAnalysis.cpp

  • issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp

    r18179 r18492  
    101101                case GiaAnalysisEnum : return new GiaAnalysis();
    102102                #endif
     103                #ifdef _HAVE_SEAICE_
     104                case SeaiceAnalysisEnum : return new SeaiceAnalysis();
     105                #endif
    103106                #ifdef _HAVE_MESHDEFORMATION_
    104107                case MeshdeformationAnalysisEnum : return new MeshdeformationAnalysis();

  • issm/trunk-jpl/src/c/analyses/analyses.h

    r18179 r18492  
    3434#include "./SmoothedSurfaceSlopeXAnalysis.h"
    3535#include "./SmoothedSurfaceSlopeYAnalysis.h"
     36#include "./SeaiceAnalysis.h"
    3637#include "./StressbalanceAnalysis.h"
    3738#include "./StressbalanceSIAAnalysis.h"

  • issm/trunk-jpl/src/c/bamg/Mesh.cpp

    r18488 r18492  
    242242                        else if (Gh.NbRef==0) delete &Gh;
    243243                }
    244                 if(&BTh && (&BTh != this)){
     244                if(&BTh && (&BTh != this) && false){
    245245                        if (BTh.NbRef>0) BTh.NbRef--;
    246246                        else if (BTh.NbRef==0) delete &BTh;
    … …  
    45054505                int verbose=bamgopts->verbose;
    45064506
    4507                 Gh.NbRef++;// add a ref to Gh
    4508 
    4509                 /*************************************************************************
    4510                  * method in 2 steps
    4511                  * 1 - compute the number of new edges to allocate
    4512                  * 2 - construct the edges
    4513                  * remark:
    4514                  * in this part we suppose to have a background mesh with the same geometry
    4515                  *
    4516                  * To construct the discretization of the new mesh we have to
    4517                  * rediscretize the boundary of background Mesh
    4518                  * because we have only the pointeur from the background mesh to the geometry.
    4519                  * We need the abcisse of the background mesh vertices on geometry
    4520                  * so a vertex is
    4521                  * 0 on GeomVertex ;
    4522                  * 1 on GeomEdge + abcisse
    4523                  * 2 internal
    4524                  *************************************************************************/
    4525 
    4526                 //Check that background mesh and current mesh do have the same geometry
    4527                 _assert_(&BTh.Gh==&Gh);
    4528                 BTh.NbRef++; // add a ref to BackGround Mesh
    4529 
    4530                 //Initialize new mesh
    4531                 BTh.SetVertexFieldOn();
    4532                 int* bcurve = new int[Gh.nbcurves]; //
    4533 
    4534                 /* There are 2 ways to make the loop
    4535                  * 1) on the geometry
    4536                  * 2) on the background mesh
    4537                  *  if you do the loop on geometry, we don't have the pointeur on background,
    4538                  *  and if you do the loop in background we have the pointeur on geometry
    4539                  * so do the walk on  background */
    4540 
    4541                 NbVerticesOnGeomVertex=0;
    4542                 NbVerticesOnGeomEdge=0;
    4543 
    4544                 /*STEP 1 copy of Required vertices*/
    4545 
    4546                 int i;
    4547                 for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
    4548                 printf("\n");
    4549                 if(NbVerticesOnGeomVertex >= maxnbv){
    4550                         _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
    4551                 }
    4552 
    4553                 VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
    4554                 VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
    4555 
    4556                 //At this point there is NO vertex but vertices should have been allocated by Init
    4557                 _assert_(vertices);
     4507                /*Intermediaries*/
     4508                int                i,k;
     4509                int                nbcurves    = 0;
     4510                int                NbNewPoints,NbEdgeCurve;
     4511                double             lcurve,lstep,s;
     4512                const int          MaxSubEdge  = 10;
     4513
     4514                R2          AB;
     4515                GeomVertex *a, *b;
     4516                BamgVertex *va,*vb;
     4517                GeomEdge   *e;
     4518
     4519                // add a ref to GH to make sure that it is not destroyed by mistake
     4520                Gh.NbRef++;
     4521
     4522                //build background mesh flag (1 if background, else 0)
     4523                bool background=(&BTh != this);
     4524
     4525                /*Build VerticesOnGeomVertex*/
     4526
     4527                //Compute the number of geometrical vertices that we are going to use to mesh
    45584528                for (i=0;i<Gh.nbv;i++){
    4559                         if (Gh[i].Required()) {//Gh vertices Required
    4560                                 vertices[nbv]  =Gh[i];
    4561                                 vertices[nbv].i=I2(0,0);
    4562                                 Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th
    4563                                 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
     4529                        if (Gh[i].Required()) NbVerticesOnGeomVertex++;
     4530                }
     4531                //allocate
     4532                VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 
     4533                if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
     4534                _assert_(nbv==0);
     4535                //Build VerticesOnGeomVertex
     4536                for (i=0;i<Gh.nbv;i++){
     4537                        /* Add vertex only if required*/
     4538                        if (Gh[i].Required()) {//Gh  vertices Required
     4539
     4540                                //Add the vertex
     4541                                _assert_(nbv<maxnbv);
     4542                                vertices[nbv]=Gh[i];
     4543
     4544                                //Add pointer from geometry (Gh) to vertex from mesh (Th)
     4545                                Gh[i].MeshVertexHook=vertices+nbv;
     4546
     4547                                //Build VerticesOnGeomVertex for current point
     4548                                VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
     4549
     4550                                //nbv increment
    45644551                                nbv++;
    45654552                        }
    4566                         else Gh[i].MeshVertexHook=0;
    4567                 }
    4568                 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){
    4569                         VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
    4570                         if (vog.IsRequiredVertex()){
    4571                                 GeomVertex* gv=vog;
    4572                                 BamgVertex *bv = vog;
    4573                                 _assert_(gv->MeshVertexHook); // use of Geom -> Th
    4574                                 VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv);
    4575                                 gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh
    4576                         }
    4577                 }
    4578                 _assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /*This might be due to MaxCornerAngle too small*/
    4579 
    4580                 /*STEP 2: reseed boundary edges*/
    4581 
    4582                 //  find the begining of the curve in BTh
    4583                 Gh.UnMarkEdges();       
    4584                 int bfind=0;
    4585                 for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1;
    4586 
    4587                 /*Loop over the backgrounf mesh BTh edges*/
    4588                 for (int iedge=0;iedge<BTh.nbe;iedge++){     
    4589                         Edge &ei = BTh.edges[iedge];
    4590 
    4591                         /*Loop over the 2 vertices of the current edge*/
    4592                         for(int je=0;je<2;je++){
    4593 
    4594                                 /* If one of the vertex is required we are in a new curve*/
    4595                                 if (ei[je].GeomEdgeHook->IsRequiredVertex()){
    4596 
    4597                                         /*Get curve number*/
    4598                                         int nc=ei.GeomEdgeHook->CurveNumber;
    4599 
    4600                                         //_printf_("Dealing with curve number " << nc << "\n");
    4601                                         //_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n");
    4602                                         //if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
    4603                                         //      _printf_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no\n");
    4604                                         //      _printf_("Do we have the right extremity? curve last  vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no\n");
    4605                                         //}
    4606                                         //BUG FIX from original bamg
    4607                                         /*Check that we are on the same edge and right vertex (0 or 1) */
    4608                                         if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge  && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
    4609                                                 bcurve[nc]=iedge*2+je;
    4610                                                 bfind++;       
    4611                                         }
    4612                                         else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge  && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
    4613                                                 bcurve[nc]=iedge*2+je;
    4614                                                 bfind++;       
    4615                                         }
    4616                                 }
    4617                         }
    4618                 }
    4619                 if (bfind!=Gh.nbcurves){
    4620                         delete [] bcurve;
    4621                         _error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)");
    4622                 }
    4623 
    4624                 // method in 2 + 1 step
    4625                 //  0.0) compute the length and the number of vertex to do allocation
    4626                 //  1.0) recompute the length
    4627                 //  1.1) compute the  vertex
    4628 
    4629                 long nbex=0,NbVerticesOnGeomEdgex=0;
    4630                 for (int step=0; step <2;step++){
    4631 
    4632                         long NbOfNewPoints=0;
    4633                         long NbOfNewEdge=0;
    4634                         long iedge;
     4553                }
     4554
     4555                /*Build VerticesOnGeomEdge*/
     4556
     4557                //check that edges is still empty (Init)
     4558                _assert_(!edges);
     4559
     4560                /* Now we are going to create the first edges corresponding
     4561                 * to the one present in the geometry provided.
     4562                 * We proceed in 2 steps
     4563                 *  -step 0: we count all the edges
     4564                 *           we allocate the number of edges at the end of step 0
     4565                 *  -step 1: the edges are created */
     4566                for (int step=0;step<2;step++){
     4567
     4568                        //initialize number of edges and number of edges max
     4569                        long nbex=0;
     4570                        nbe=0;
     4571                        long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
    46354572                        Gh.UnMarkEdges();       
    4636                         double L=0;
    4637 
    4638                         /*Go through all geometrical curve*/
    4639                         for (int icurve=0;icurve<Gh.nbcurves;icurve++){
    4640 
    4641                                 /*Get edge and vertex (index) of background mesh on this curve*/
    4642                                 iedge=bcurve[icurve]/2;
    4643                                 int jedge=bcurve[icurve]%2;
    4644 
    4645                                 /*Get edge of Bth with index iedge*/
    4646                                 Edge &ei = BTh.edges[iedge];
    4647 
    4648                                 /*Initialize variables*/
    4649                                 double Lstep=0;             // step between two points   (phase==1)
    4650                                 long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
    4651 
    4652                                 /*Do phase 0 to step*/
    4653                                 for(int phase=0;phase<=step;phase++){
    4654 
    4655                                         /*Current curve pointer*/
    4656                                         Curve *curve= Gh.curves+icurve;
    4657 
    4658                                         /*Get index of current curve*/
    4659                                         int icurveequi= Gh.GetId(curve);
    4660 
    4661                                         /*For phase 0, check that we are at the begining of the curve only*/
    4662                                         if(phase==0 &&  icurveequi!=icurve)  continue;
    4663 
    4664                                         int   k0=jedge,k1;
    4665                                         Edge* pe=  BTh.edges+iedge;
    4666                                         int   iedgeequi=bcurve[icurveequi]/2;
    4667                                         int   jedgeequi=bcurve[icurveequi]%2;
    4668 
    4669                                         int k0equi=jedgeequi,k1equi;             
    4670                                         Edge * peequi= BTh.edges+iedgeequi;
    4671                                         GeomEdge *ongequi = peequi->GeomEdgeHook;
    4672 
    4673                                         double sNew=Lstep;// abscisse of the new points (phase==1)
    4674                                         L=0;// length of the curve
    4675                                         long i=0;// index of new points on the curve
    4676                                         GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook;
    4677                                         BamgVertex *A0;
    4678                                         A0 = GA0->MeshVertexHook;  // the vertex in new mesh
    4679                                         BamgVertex *A1;
    4680                                         VertexOnGeom *GA1;
    4681                                         Edge* PreviousNewEdge = 0;
    4682 
    4683                                         // New Curve phase
    4684                                         _assert_(A0-vertices>=0 && A0-vertices<nbv);
    4685                                         if(ongequi->Required()){
    4686                                                 GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook;
    4687                                                 A1 = GA1->MeshVertexHook;  //
    4688                                         }       
    4689                                         else {
    4690                                                 for(;;){
    4691                                                         Edge &ee=*pe;
    4692                                                         Edge &eeequi=*peequi;
    4693                                                         k1 = 1-k0; // next vertex of the edge
    4694                                                         k1equi= 1 - k0equi;
    4695                                                         _assert_(pe && ee.GeomEdgeHook);
    4696                                                         ee.GeomEdgeHook->SetMark();
    4697                                                         BamgVertex & v0=ee[0], & v1=ee[1];
    4698                                                         R2 AB=(R2)v1-(R2)v0;
    4699                                                         double L0=L,LAB;
    4700                                                         LAB=LengthInterpole(v0.m,v1.m,AB);
    4701                                                         L+= LAB;
    4702 
    4703                                                         if (phase){
    4704                                                                 // computation of the new points for the given curve
    4705                                                                 while ((i!=NbCreatePointOnCurve) && sNew<=L) {
    4706 
    4707                                                                         //some checks
    4708                                                                         _assert_(sNew>=L0);
    4709                                                                         _assert_(LAB);
    4710                                                                         _assert_(vertices && nbv<maxnbv);
    4711                                                                         _assert_(edges && nbe<nbex);
    4712                                                                         _assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
    4713 
    4714                                                                         // new vertex on edge
    4715                                                                         A1=vertices+nbv++;
    4716                                                                         GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
    4717                                                                         Edge* e = edges + nbe++;
    4718                                                                         double se= (sNew-L0)/LAB;
    4719                                                                         if (se<0 || se>=1.000000001){
    4720                                                                                 _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
    4721                                                                         }
    4722                                                                         se = abscisseInterpole(v0.m,v1.m,AB,se,1);
    4723                                                                         if (se<0 || se>1){
    4724                                                                                 _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
    4725                                                                         }
    4726                                                                         se = k1         ? se : 1. - se;
    4727                                                                         se = k1==k1equi ? se : 1. - se;
    4728                                                                         VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
    4729                                                                         ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
    4730                                                                         A1->ReferenceNumber = eeequi.ReferenceNumber;
    4731                                                                         A1->DirOfSearch =NoDirOfSearch;
    4732                                                                         e->GeomEdgeHook = ongequi;
    4733                                                                         e->v[0]=A0;
    4734                                                                         e->v[1]=A1;
    4735                                                                         e->ReferenceNumber = eeequi.ReferenceNumber;
    4736                                                                         e->adj[0]=PreviousNewEdge;
    4737 
    4738                                                                         if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
    4739                                                                         PreviousNewEdge=e;
    4740                                                                         A0=A1;
    4741                                                                         sNew += Lstep;
    4742                                                                         if (++i== NbCreatePointOnCurve) break;
     4573                        nbcurves=0;
     4574
     4575                        //go through the edges of the geometry
     4576                        for (i=0;i<Gh.nbe;i++){
     4577
     4578                                //ei = current Geometrical edge
     4579                                GeomEdge &ei=Gh.edges[i];   
     4580
     4581                                //loop over the two vertices of the edge ei
     4582                                for(int j=0;j<2;j++) {
     4583
     4584                                        /*Take only required vertices (corner->beginning of a new curve)*/
     4585                                        if (!ei.Mark() && ei[j].Required()){
     4586
     4587                                                long  nbvend=0;
     4588                                                Edge* PreviousNewEdge=NULL;
     4589                                                lstep = -1;
     4590
     4591                                                /*If Edge is required (do that only once for the 2 vertices)*/
     4592                                                if(ei.Required()){
     4593                                                        if (j==0){
     4594                                                                //do not create internal points if required (take it as is)
     4595                                                                if(step==0) nbe++;
     4596                                                                else{
     4597                                                                        e=&ei;
     4598                                                                        a=ei(0);
     4599                                                                        b=ei(1);
     4600
     4601                                                                        //check that edges has been allocated
     4602                                                                        _assert_(edges);
     4603                                                                        edges[nbe].v[0]=a->MeshVertexHook;
     4604                                                                        edges[nbe].v[1]=b->MeshVertexHook;;
     4605                                                                        edges[nbe].ReferenceNumber = e->ReferenceNumber;
     4606                                                                        edges[nbe].GeomEdgeHook = e;
     4607                                                                        edges[nbe].adj[0] = 0;
     4608                                                                        edges[nbe].adj[1] = 0;
     4609                                                                        nbe++;
    47434610                                                                }
    47444611                                                        }
    4745 
    4746                                                         //some checks
    4747                                                         _assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
    4748                                                         if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
    4749                                                                 _assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
    4750                                                                 GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
    4751                                                                 A1=GA1->MeshVertexHook;// the vertex in new mesh
    4752                                                                 _assert_(A1-vertices>=0 && A1-vertices<nbv);
    4753                                                                 break;
     4612                                                }
     4613
     4614                                                /*If Edge is not required: we are on a curve*/
     4615                                                else {
     4616                                                        for (int kstep=0;kstep<=step;kstep++){
     4617                                                                //kstep=0, compute number of edges (discretize curve)
     4618                                                                //kstep=1  create the points and edge
     4619                                                                PreviousNewEdge=0;
     4620                                                                NbNewPoints=0;
     4621                                                                NbEdgeCurve=0;
     4622                                                                if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
     4623                                                                lcurve =0;
     4624                                                                s = lstep; //-1 initially, then length of each sub edge
     4625
     4626                                                                /*reminder: i = edge number, j=[0;1] vertex index in edge*/
     4627                                                                k=j;            // k = vertex index in edge (0 or 1)
     4628                                                                e=&ei;          // e = reference of current edge
     4629                                                                a=ei(k);        // a = pointer toward the kth vertex of the current edge
     4630                                                                va = a->MeshVertexHook; // va = pointer toward mesh vertex associated
     4631                                                                e->SetMark();   // Mark edge
     4632
     4633                                                                /*Loop until we reach the end of the curve*/
     4634                                                                for(;;){
     4635                                                                        k = 1-k;            // other vertx index of the curve
     4636                                                                        b = (*e)(k);        // b = pointer toward the other vertex of the current edge
     4637                                                                        AB= b->r - a->r;   // AB = vector of the current edge
     4638                                                                        Metric MA = background ? BTh.MetricAt(a->r) :a->m ;  //Get metric associated to A
     4639                                                                        Metric MB = background ? BTh.MetricAt(b->r) :b->m ;  //Get metric associated to B
     4640                                                                        double ledge = (MA(AB) + MB(AB))/2;                  //Get edge length in metric
     4641
     4642                                                                        /* We are now creating the mesh edges from the geometrical edge selected above.
     4643                                                                         * The edge will be divided according to the metric previously computed and cannot
     4644                                                                         * be divided more than 10 times (MaxSubEdge). */
     4645
     4646                                                                        //By default, there is only one subedge that is the geometrical edge itself
     4647                                                                        int NbSubEdge = 1;
     4648
     4649                                                                        //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
     4650                                                                        double lSubEdge[MaxSubEdge];
     4651
     4652                                                                        //Build Subedges according to the edge length
     4653                                                                        if (ledge < 1.5){
     4654                                                                                //if ledge < 1.5 (between one and 2), take the edge as is
     4655                                                                                lSubEdge[0] = ledge;
     4656                                                                        }
     4657                                                                        //else, divide the edge
     4658                                                                        else {
     4659                                                                                //compute number of subedges (division of the edge), Maximum is 10
     4660                                                                                NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
     4661                                                                                /*Now, we are going to divide the edge according to the metric.
     4662                                                                                 * Get segment by sement along the edge.
     4663                                                                                 * Build lSubEdge, which holds the distance between the first vertex
     4664                                                                                 * of the edge and the next point on the edge according to the
     4665                                                                                 * discretization (each SubEdge is AB)*/
     4666                                                                                R2 A,B;
     4667                                                                                A=a->r;
     4668                                                                                Metric MAs=MA,MBs;
     4669                                                                                ledge=0;
     4670                                                                                double x =0, xstep= 1./NbSubEdge;
     4671                                                                                for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
     4672                                                                                        x += xstep;
     4673                                                                                        B =  e->F(k ? x : 1-x);
     4674                                                                                        MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
     4675                                                                                        AB = A-B;
     4676                                                                                        lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2);
     4677                                                                                }
     4678                                                                        }
     4679
     4680                                                                        double lcurveb = lcurve+ledge;
     4681
     4682                                                                        /*Now, create corresponding points*/
     4683                                                                        while(s>=lcurve && s<=lcurveb && nbv<nbvend){
     4684
     4685                                                                                /*Schematic of current curve
     4686                                                                                 *
     4687                                                                                 *  a                   vb                  b          // vertex
     4688                                                                                 *  0              ll0     ll1              ledge      // length from a
     4689                                                                                 *  + --- + - ... - + --S-- + --- + - ... - +          // where is S
     4690                                                                                 *  0              kk0     kk1              NbSubEdge  // Sub edge index
     4691                                                                                 *
     4692                                                                                 */
     4693
     4694                                                                                double ss = s-lcurve;
     4695
     4696                                                                                /*Find the SubEdge containing ss using Dichotomy*/
     4697                                                                                int kk0=-1,kk1=NbSubEdge-1,kkk;
     4698                                                                                double ll0=0,ll1=ledge,llk;
     4699                                                                                while (kk1-kk0>1){
     4700                                                                                        if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
     4701                                                                                         kk1=kkk,ll1=llk;
     4702                                                                                        else
     4703                                                                                         kk0=kkk,ll0=llk;
     4704                                                                                }
     4705                                                                                _assert_(kk1!=kk0);
     4706
     4707                                                                                /*Curvilinear coordinate in [0 1] of ss in current edge*/
     4708                                                                                // WARNING: This is what we would do
     4709                                                                                // ssa = (ss-ll0)/(ll1-ll0);
     4710                                                                                // aa = (kk0+ssa)/NbSubEdge
     4711                                                                                // This is what Bamg does:
     4712                                                                                double sbb = (ss-ll0)/(ll1-ll0);
     4713                                                                                /*Curvilinear coordinate in [0 1] of ss in current curve*/
     4714                                                                                double bb = (kk1+sbb)/NbSubEdge;
     4715                                                                                double aa = 1-bb;
     4716
     4717                                                                                // new vertex on edge
     4718                                                                                vb = &vertices[nbv++];
     4719                                                                                vb->m = Metric(aa,a->m,bb,b->m);
     4720                                                                                vb->ReferenceNumber = e->ReferenceNumber;
     4721                                                                                vb->DirOfSearch =NoDirOfSearch;
     4722                                                                                double abcisse = k ? bb : aa;
     4723                                                                                vb->r =  e->F(abcisse);
     4724                                                                                VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);       
     4725
     4726                                                                                // to take into account the direction of the edge
     4727                                                                                s += lstep;
     4728                                                                                edges[nbe].v[0]=va;
     4729                                                                                edges[nbe].v[1]=vb;
     4730                                                                                edges[nbe].ReferenceNumber =e->ReferenceNumber;
     4731                                                                                edges[nbe].GeomEdgeHook = e;
     4732                                                                                edges[nbe].adj[0] = PreviousNewEdge;
     4733                                                                                if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
     4734                                                                                PreviousNewEdge=edges+nbe;
     4735                                                                                nbe++;
     4736                                                                                va = vb;
     4737                                                                        }
     4738
     4739                                                                        /*We just added one edge to the curve: Go to the next one*/
     4740                                                                        lcurve = lcurveb;
     4741                                                                        e->SetMark();
     4742                                                                        a=b;
     4743
     4744                                                                        /*If b is required, we are on a new curve->break*/
     4745                                                                        if (b->Required()) break;
     4746                                                                        int kprev=k;
     4747                                                                        k = e->AdjVertexIndex[kprev];// next vertices
     4748                                                                        e = e->Adj[kprev];
     4749                                                                        _assert_(e);
     4750                                                                }// for(;;)
     4751                                                                vb = b->MeshVertexHook;
     4752
     4753                                                                /*Number of edges in the last disretized curve*/
     4754                                                                NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
     4755                                                                /*Number of internal vertices in the last disretized curve*/
     4756                                                                NbNewPoints = NbEdgeCurve-1;
     4757                                                                if(!kstep){
     4758                                                                        NbVerticesOnGeomEdge0 += NbNewPoints;
     4759                                                                        nbcurves++;
     4760                                                                }
     4761                                                                nbvend=nbv+NbNewPoints;
     4762                                                                lstep = lcurve / NbEdgeCurve; //approximately one
     4763                                                        }// end of curve --
     4764                                                        if (edges) { // last edges of the curves
     4765                                                                edges[nbe].v[0]=va;
     4766                                                                edges[nbe].v[1]=vb;
     4767                                                                edges[nbe].ReferenceNumber = e->ReferenceNumber;
     4768                                                                edges[nbe].GeomEdgeHook = e;
     4769                                                                edges[nbe].adj[0] = PreviousNewEdge;
     4770                                                                edges[nbe].adj[1] = 0;
     4771                                                                if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];
     4772                                                                nbe++;
    47544773                                                        }
    4755                                                         if (!ee.adj[k1]) {
    4756                                                                 _error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices);
    4757                                                         }
    4758                                                         pe = ee.adj[k1]; // next edge
    4759                                                         k0 = pe->Intersection(ee);
    4760                                                         peequi= eeequi.adj[k1equi];  // next edge
    4761                                                         k0equi=peequi->Intersection(eeequi);           
    4762                                                 }// for(;;) end of the curve
     4774                                                        else nbe += NbEdgeCurve;
     4775                                                } // end on  curve ---
    47634776                                        }
    4764 
    4765                                         if (phase){ // construction of the last edge
    4766                                                 Edge* e=edges + nbe++;
    4767                                                 e->GeomEdgeHook  = ongequi;
    4768                                                 e->v[0]=A0;
    4769                                                 e->v[1]=A1;
    4770                                                 e->ReferenceNumber = peequi->ReferenceNumber;
    4771                                                 e->adj[0]=PreviousNewEdge;
    4772                                                 e->adj[1]=0;
    4773                                                 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
    4774                                                 PreviousNewEdge = e;
    4775 
    4776                                                 _assert_(i==NbCreatePointOnCurve);
    4777                                         }
    4778 
    4779                                         if (!phase)  { //
    4780                                                 long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
    4781                                                 Lstep = L/NbSegOnCurve;
    4782                                                 NbCreatePointOnCurve = NbSegOnCurve-1;
    4783                                                 NbOfNewEdge += NbSegOnCurve;
    4784                                                 NbOfNewPoints += NbCreatePointOnCurve;
    4785                                         }
    4786                                 }
    4787                         }//  end of curve loop
    4788 
    4789                         //Allocate memory
    4790                         if(step==0){
    4791                                 if(nbv+NbOfNewPoints > maxnbv) {
    4792                                         _error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv);
    4793                                 }
    4794                                 edges = new Edge[NbOfNewEdge];
    4795                                 nbex = NbOfNewEdge;
    4796                                 if(NbOfNewPoints) {
    4797                                         VerticesOnGeomEdge    = new VertexOnGeom[NbOfNewPoints];
    4798                                         NbVertexOnBThEdge     = NbOfNewPoints;
    4799                                         VertexOnBThEdge       = new  VertexOnEdge[NbOfNewPoints];
    4800                                         NbVerticesOnGeomEdgex = NbOfNewPoints;
    4801                                 }
    4802                                 NbOfNewPoints =0;
    4803                                 NbOfNewEdge = 0;
    4804                         }
    4805                 }
    4806                 _assert_(nbe!=0);
    4807                 delete [] bcurve;
     4777                                }
     4778                        } // for (i=0;i<nbe;i++)
     4779                        if(!step) {
     4780                                _assert_(!edges);
     4781                                _assert_(!VerticesOnGeomEdge);
     4782
     4783                                edges = new Edge[nbex=nbe];
     4784                                if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
     4785
     4786                                // do the vertex on a geometrical vertex
     4787                                _assert_(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0);
     4788                                NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;       
     4789                        }
     4790                        else{
     4791                                _assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
     4792                        }
     4793                }
     4794
    48084795
    48094796                //Insert points inside existing triangles

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r18237 r18492  
    2020
    2121        bool isefficientlayer;
    22         int  hydrology_model,smb_model;
     22        int  hydrology_model,smb_model,materials_type;
    2323        iomodel->Constant(&hydrology_model,HydrologyModelEnum);
    2424        iomodel->Constant(&smb_model,SurfaceforcingsEnum);
     25        iomodel->Constant(&materials_type,MaterialsEnum);
    2526
    2627        this->mid = matpar_mid;
    27         iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
    28         iomodel->Constant(&this->rho_water,MaterialsRhoSeawaterEnum);
    29         iomodel->Constant(&this->rho_freshwater,MaterialsRhoFreshwaterEnum);
    30         iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum);
    31         iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
    32         iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum);
    33         iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum);
    34         iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum);
    35         iomodel->Constant(&this->beta,MaterialsBetaEnum);
    36         iomodel->Constant(&this->meltingpoint,MaterialsMeltingpointEnum);
    37         iomodel->Constant(&this->referencetemperature,ConstantsReferencetemperatureEnum);
    38         iomodel->Constant(&this->mixed_layer_capacity,MaterialsMixedLayerCapacityEnum);
    39         iomodel->Constant(&this->thermal_exchange_velocity,MaterialsThermalExchangeVelocityEnum);
    40         iomodel->Constant(&this->g,ConstantsGEnum);
    41 
    42         switch(smb_model){
    43                 case SMBEnum:
    44                         /*Nothing to add*/
    45                         break;
    46                 case SMBpddEnum:
    47                         iomodel->Constant(&this->desfac,SurfaceforcingsDesfacEnum);
    48                         iomodel->Constant(&this->s0p,SurfaceforcingsS0pEnum);
    49                         break;
    50                 case SMBgradientsEnum:
    51                         /*Nothing to add*/
    52                         break;
    53                 case SMBhenningEnum:
    54                         /*Nothing to add*/
    55                         break;
    56                 case SMBcomponentsEnum:
    57                         /*Nothing to add*/
    58                         break;
    59                 case SMBmeltcomponentsEnum:
    60                         /*Nothing to add*/
     28
     29        switch(materials_type){
     30                case MaticeEnum:
     31                        iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
     32                        iomodel->Constant(&this->rho_water,MaterialsRhoSeawaterEnum);
     33                        iomodel->Constant(&this->rho_freshwater,MaterialsRhoFreshwaterEnum);
     34                        iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum);
     35                        iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
     36                        iomodel->Constant(&this->thermalconductivity,MaterialsThermalconductivityEnum);
     37                        iomodel->Constant(&this->temperateiceconductivity,MaterialsTemperateiceconductivityEnum);
     38                        iomodel->Constant(&this->latentheat,MaterialsLatentheatEnum);
     39                        iomodel->Constant(&this->beta,MaterialsBetaEnum);
     40                        iomodel->Constant(&this->meltingpoint,MaterialsMeltingpointEnum);
     41                        iomodel->Constant(&this->referencetemperature,ConstantsReferencetemperatureEnum);
     42                        iomodel->Constant(&this->mixed_layer_capacity,MaterialsMixedLayerCapacityEnum);
     43                        iomodel->Constant(&this->thermal_exchange_velocity,MaterialsThermalExchangeVelocityEnum);
     44                        iomodel->Constant(&this->g,ConstantsGEnum);
     45
     46                        switch(smb_model){
     47                                case SMBEnum:
     48                                        /*Nothing to add*/
     49                                        break;
     50                                case SMBpddEnum:
     51                                        iomodel->Constant(&this->desfac,SurfaceforcingsDesfacEnum);
     52                                        iomodel->Constant(&this->s0p,SurfaceforcingsS0pEnum);
     53                                        break;
     54                                case SMBgradientsEnum:
     55                                        /*Nothing to add*/
     56                                        break;
     57                                case SMBhenningEnum:
     58                                        /*Nothing to add*/
     59                                        break;
     60                                case SMBcomponentsEnum:
     61                                        /*Nothing to add*/
     62                                        break;
     63                                case SMBmeltcomponentsEnum:
     64                                        /*Nothing to add*/
     65                                        break;
     66                                default:
     67                                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
     68                        }
     69                        if(hydrology_model==HydrologydcEnum){
     70                                iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
     71                                iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);
     72                                iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum);
     73                                iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
     74                                iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     75
     76                                if(isefficientlayer){
     77                                        iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
     78                                        iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
     79                                        iomodel->Constant(&this->epl_init_thickness,HydrologydcEplInitialThicknessEnum);
     80                                        iomodel->Constant(&this->epl_max_thickness,HydrologydcEplMaxThicknessEnum);
     81                                        iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
     82                                }
     83                        }
     84                        else if(hydrology_model==HydrologyshreveEnum){
     85                                /*Nothing to add*/
     86                        }
     87                        else{
     88                                _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
     89                        }
     90
     91                        /*gia: */
     92                        iomodel->Constant(&this->lithosphere_shear_modulus,MaterialsLithosphereShearModulusEnum);
     93                        iomodel->Constant(&this->lithosphere_density,MaterialsLithosphereDensityEnum);
     94                        iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
     95                        iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);
     96                        break;
     97                case MatseaiceEnum:
     98                        iomodel->Constant(&this->poisson,MaterialsPoissonEnum);
     99                        iomodel->Constant(&this->young_modulus,MaterialsYoungModulusEnum);
     100                        iomodel->Constant(&this->ridging_exponent,MaterialsRidgingExponentEnum);
     101                        iomodel->Constant(&this->cohesion,MaterialsCohesionEnum);
     102                        iomodel->Constant(&this->internal_friction_coef,MaterialsInternalFrictionCoefEnum);
     103                        iomodel->Constant(&this->compression_coef,MaterialsCompressionCoefEnum);
     104                        iomodel->Constant(&this->traction_coef,MaterialsTractionCoefEnum);
     105                        iomodel->Constant(&this->g,ConstantsGEnum);
    61106                        break;
    62107                default:
    63                         _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    64         }
    65         if(hydrology_model==HydrologydcEnum){
    66                 iomodel->Constant(&this->sediment_compressibility,HydrologydcSedimentCompressibilityEnum);
    67                 iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);
    68                 iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum);
    69                 iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
    70                 iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    71 
    72                 if(isefficientlayer){
    73                                 iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
    74                                 iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
    75                                 iomodel->Constant(&this->epl_init_thickness,HydrologydcEplInitialThicknessEnum);
    76                                 iomodel->Constant(&this->epl_max_thickness,HydrologydcEplMaxThicknessEnum);
    77                                 iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
    78                 }
    79         }
    80         else if(hydrology_model==HydrologyshreveEnum){
    81                 /*Nothing to add*/
    82         }
    83         else{
    84                 _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");
    85         }
    86 
    87         /*gia: */
    88         iomodel->Constant(&this->lithosphere_shear_modulus,MaterialsLithosphereShearModulusEnum);
    89         iomodel->Constant(&this->lithosphere_density,MaterialsLithosphereDensityEnum);
    90         iomodel->Constant(&this->mantle_shear_modulus,MaterialsMantleShearModulusEnum);
    91         iomodel->Constant(&this->mantle_density,MaterialsMantleDensityEnum);
     108                        _error_("Material not supported yet");
     109        }
    92110}
    93111/*}}}*/
    … …  
    293311                case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
    294312                case ConstantsGEnum:                         return this->g;
     313                case MaterialsPoissonEnum:                   return this->poisson;
     314                case MaterialsYoungModulusEnum:              return this->young_modulus;
     315                case MaterialsRidgingExponentEnum:                              return this->ridging_exponent;
     316                case MaterialsCohesionEnum:                         return this->cohesion;
     317                case MaterialsInternalFrictionCoefEnum:         return this->internal_friction_coef;
     318                case MaterialsCompressionCoefEnum:                              return this->compression_coef;
     319                case MaterialsTractionCoefEnum:                                 return this->traction_coef;
    295320                default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
    296321        }

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r18237 r18492  
    5050                IssmDouble mantle_shear_modulus;
    5151                IssmDouble mantle_density;
     52
     53                /*Sea ice*/
     54                IssmDouble poisson;
     55                IssmDouble young_modulus;
     56                IssmDouble ridging_exponent;
     57                IssmDouble cohesion;
     58                IssmDouble internal_friction_coef;
     59                IssmDouble compression_coef;
     60                IssmDouble traction_coef;
    5261
    5362        public:

  • issm/trunk-jpl/src/c/classes/Vertex.cpp

    r18063 r18492  
    3232        this->domaintype     = iomodel->domaintype;
    3333
    34         _assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum) && iomodel->numbernodetoelementconnectivity);
    3534        switch(iomodel->domaintype){
    3635                case Domain3DEnum:
    37                 case Domain2DhorizontalEnum:
     36                        _assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum));
    3837                        this->sigma = (iomodel->Data(MeshZEnum)[i]-iomodel->Data(BaseEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
    3938                        break;
     39                case Domain2DhorizontalEnum:
     40                        this->sigma = 0.;
     41                        break;
    4042                case Domain2DverticalEnum:
     43                        _assert_(iomodel->Data(BaseEnum) && iomodel->Data(ThicknessEnum));
    4144                        this->sigma = (iomodel->Data(MeshYEnum)[i]-iomodel->Data(BaseEnum)[i])/(iomodel->Data(ThicknessEnum)[i]);
    4245                        break;
    4346        }
    4447
     48        _assert_(iomodel->numbernodetoelementconnectivity);
    4549        this->connectivity = iomodel->numbernodetoelementconnectivity[i];
    4650

  • issm/trunk-jpl/src/c/classes/classes.h

    r17472 r18492  
    7878#include "./Materials/Material.h"
    7979#include "./Materials/Matice.h"
     80#include "./Materials/Matseaice.h"
    8081#include "./Materials/Matpar.h"
    8182

  • issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp

    r18485 r18492  
    148148                        break;
    149149
     150                case SeaiceSolutionEnum:
     151                        numanalyses=1;
     152                        analyses=xNew<int>(numanalyses);
     153                        analyses[0]=SeaiceAnalysisEnum;
     154                        break;
     155
    150156                default:
    151157                        _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");

  • issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp

    r17831 r18492  
    6969                        solutioncore=&damage_core;
    7070                        break;
     71                case SeaiceSolutionEnum:
     72                        solutioncore=&seaice_core;
     73                        break;
    7174                default:
    7275                        _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");

  • issm/trunk-jpl/src/c/cores/cores.h

    r18134 r18492  
    4545void gia_core(FemModel* femmodel);
    4646void damage_core(FemModel* femmodel);
     47void seaice_core(FemModel* femmodel);
    4748IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
    4849

  • issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp

    r15104 r18492  
    200200
    201201                /*clean up*/
     202                printf("-------------- file: Bamgx.cpp line: %i\n",__LINE__);
    202203                delete &Th;
     204                printf("-------------- file: Bamgx.cpp line: %i\n",__LINE__);
    203205                //delete &BTh;
    204206                /*}}}*/
    … …  
    207209        /*No error return*/
    208210        if (verbosity>1) _printf_("   Exiting Bamg.\n");
     211        printf("-------------- file: Bamgx.cpp line: %i\n",__LINE__);
    209212        return noerr;
    210213

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r18237 r18492  
    7979                        }
    8080                        break;
     81                case MatseaiceEnum:
     82                        iomodel->FetchDataToInput(elements,MaterialsDamageEnum);
     83                        for (i=0;i<iomodel->numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matseaice(i+1,i,iomodel));
     84                        break;
    8185                default:
    8286                        _error_("Materials "<<EnumToStringx(materials_type)<<" not supported");

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r18093 r18492  
    4343        parameters->AddObject(iomodel->CopyConstantObject(DomainTypeEnum));
    4444        parameters->AddObject(iomodel->CopyConstantObject(DomainDimensionEnum));
    45         parameters->AddObject(iomodel->CopyConstantObject(MeshElementtypeEnum));
    4645        parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum));
    47         parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
    48         parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum));
    4946        parameters->AddObject(iomodel->CopyConstantObject(ConstantsYtsEnum));
    5047        parameters->AddObject(iomodel->CopyConstantObject(TimesteppingStartTimeEnum));
    … …  
    5754        parameters->AddObject(iomodel->CopyConstantObject(DebugProfilingEnum));
    5855        parameters->AddObject(iomodel->CopyConstantObject(MeshAverageVertexConnectivityEnum));
    59         parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
    6056        parameters->AddObject(iomodel->CopyConstantObject(SettingsWaitonlockEnum));
    6157        parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofelementsEnum));
    … …  
    6359        parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsOnNodesEnum));
    6460        parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum));
    65         parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
    66         parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum));
    67         parameters->AddObject(iomodel->CopyConstantObject(TransientIsmasstransportEnum));
    68         parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
    69         parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
    70         parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum));
    71         parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
    72         parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
    73         parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
    74         parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
    7561        parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));
    7662        parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
    7763        parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
    7864        parameters->AddObject(iomodel->CopyConstantObject(InversionTypeEnum));
    79         parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
    80 
    81         /*For stress balance only*/
    82         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
    83         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
    84         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
    85         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
    86         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
    87         parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
    88         if(iomodel->domaintype==Domain3DEnum)
    89          parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));
     65        if(solution_type==SeaiceSolutionEnum){
     66        }
     67        else{
     68                parameters->AddObject(iomodel->CopyConstantObject(MeshElementtypeEnum));
     69                parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
     70                parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum));
     71                parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
     72                parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
     73                parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum));
     74                parameters->AddObject(iomodel->CopyConstantObject(TransientIsmasstransportEnum));
     75                parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
     76                parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
     77                parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum));
     78                parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum));
     79                parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
     80                parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
     81                parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
     82                parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));
     83
     84                /*For stress balance only*/
     85                parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
     86                parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
     87                parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
     88                parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
     89                parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
     90                parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
     91                if(iomodel->domaintype==Domain3DEnum)
     92                 parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));
     93        }
    9094
    9195        /*Surface mass balance parameters*/
    … …  
    129133                        /*Nothing to add to parameters*/
    130134                        break;
     135                case SeaiceatmEnum:
     136                        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsRhoAirEnum));
     137                        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsAirCoefEnum));
     138                        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsAirLinDragCoefEnum));
     139                        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsAirQuadDragCoefEnum));
     140                        break;
    131141                default:
    132142                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    … …  
    144154                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsDeepwaterElevationEnum));
    145155                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsUpperwaterElevationEnum));
     156                        break;
     157                case SeaiceoceanEnum:
     158                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsRhoOceanEnum));
     159                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanCoefEnum));
     160                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanLinDragCoefEnum));
     161                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanQuadDragCoefEnum));
     162                        parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsOceanTurningAngleEnum));
    146163                        break;
    147164                default:

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18476 r18492  
    387387        GiaSolutionEnum,
    388388        GiaAnalysisEnum,
     389        SeaiceSolutionEnum,
     390        SeaiceAnalysisEnum,
    389391        MeshdeformationSolutionEnum,
    390392        MeshdeformationAnalysisEnum,
    … …  
    742744        LevelsetfunctionPicardEnum,
    743745        /*}}}*/
     746        /*Sea Ice {{{*/
     747        SeaiceatmEnum,
     748        SeaiceoceanEnum,
     749        SeaiceThicknessEnum,
     750        SeaiceConcentrationEnum,
     751        SeaiceSpcvxEnum,
     752        SeaiceSpcvyEnum,
     753        BasalforcingsRhoOceanEnum,
     754        BasalforcingsOceanCoefEnum,
     755        BasalforcingsOceanLinDragCoefEnum,
     756        BasalforcingsOceanQuadDragCoefEnum,
     757        BasalforcingsOceanTurningAngleEnum,
     758        BasalforcingsOceanSshEnum,
     759        BasalforcingsOceanVxEnum,
     760        BasalforcingsOceanVyEnum,
     761        SurfaceforcingsRhoAirEnum,
     762        SurfaceforcingsAirCoefEnum,
     763        SurfaceforcingsAirLinDragCoefEnum,
     764        SurfaceforcingsAirQuadDragCoefEnum,
     765        SurfaceforcingsWindVxEnum,
     766        SurfaceforcingsWindVyEnum,
     767        MatseaiceEnum,
     768        MaterialsPoissonEnum,
     769        MaterialsYoungModulusEnum,
     770        MaterialsDamageEnum,
     771        MaterialsRidgingExponentEnum,
     772        MaterialsCohesionEnum,
     773        MaterialsInternalFrictionCoefEnum,
     774        MaterialsCompressionCoefEnum,
     775        MaterialsTractionCoefEnum,
     776        VxStarEnum,
     777        VyStarEnum,
     778        /*}}}*/
    744779        MaximumNumberOfDefinitionsEnum
    745780};

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18476 r18492  
    390390                case GiaSolutionEnum : return "GiaSolution";
    391391                case GiaAnalysisEnum : return "GiaAnalysis";
     392                case SeaiceSolutionEnum : return "SeaiceSolution";
     393                case SeaiceAnalysisEnum : return "SeaiceAnalysis";
    392394                case MeshdeformationSolutionEnum : return "MeshdeformationSolution";
    393395                case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
    … …  
    702704                case LevelsetfunctionSlopeYEnum : return "LevelsetfunctionSlopeY";
    703705                case LevelsetfunctionPicardEnum : return "LevelsetfunctionPicard";
     706                case SeaiceatmEnum : return "Seaiceatm";
     707                case SeaiceoceanEnum : return "Seaiceocean";
     708                case SeaiceThicknessEnum : return "SeaiceThickness";
     709                case SeaiceConcentrationEnum : return "SeaiceConcentration";
     710                case SeaiceSpcvxEnum : return "SeaiceSpcvx";
     711                case SeaiceSpcvyEnum : return "SeaiceSpcvy";
     712                case BasalforcingsRhoOceanEnum : return "BasalforcingsRhoOcean";
     713                case BasalforcingsOceanCoefEnum : return "BasalforcingsOceanCoef";
     714                case BasalforcingsOceanLinDragCoefEnum : return "BasalforcingsOceanLinDragCoef";
     715                case BasalforcingsOceanQuadDragCoefEnum : return "BasalforcingsOceanQuadDragCoef";
     716                case BasalforcingsOceanTurningAngleEnum : return "BasalforcingsOceanTurningAngle";
     717                case BasalforcingsOceanSshEnum : return "BasalforcingsOceanSsh";
     718                case BasalforcingsOceanVxEnum : return "BasalforcingsOceanVx";
     719                case BasalforcingsOceanVyEnum : return "BasalforcingsOceanVy";
     720                case SurfaceforcingsRhoAirEnum : return "SurfaceforcingsRhoAir";
     721                case SurfaceforcingsAirCoefEnum : return "SurfaceforcingsAirCoef";
     722                case SurfaceforcingsAirLinDragCoefEnum : return "SurfaceforcingsAirLinDragCoef";
     723                case SurfaceforcingsAirQuadDragCoefEnum : return "SurfaceforcingsAirQuadDragCoef";
     724                case SurfaceforcingsWindVxEnum : return "SurfaceforcingsWindVx";
     725                case SurfaceforcingsWindVyEnum : return "SurfaceforcingsWindVy";
     726                case MatseaiceEnum : return "Matseaice";
     727                case MaterialsPoissonEnum : return "MaterialsPoisson";
     728                case MaterialsYoungModulusEnum : return "MaterialsYoungModulus";
     729                case MaterialsDamageEnum : return "MaterialsDamage";
     730                case MaterialsRidgingExponentEnum : return "MaterialsRidgingExponent";
     731                case MaterialsCohesionEnum : return "MaterialsCohesion";
     732                case MaterialsInternalFrictionCoefEnum : return "MaterialsInternalFrictionCoef";
     733                case MaterialsCompressionCoefEnum : return "MaterialsCompressionCoef";
     734                case MaterialsTractionCoefEnum : return "MaterialsTractionCoef";
     735                case VxStarEnum : return "VxStar";
     736                case VyStarEnum : return "VyStar";
    704737                case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
    705738                default : return "unknown";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18476 r18492  
    399399              else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
    400400              else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
     401              else if (strcmp(name,"SeaiceSolution")==0) return SeaiceSolutionEnum;
     402              else if (strcmp(name,"SeaiceAnalysis")==0) return SeaiceAnalysisEnum;
    401403              else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
    402404              else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
    … …  
    504506              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    505507              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    506               else if (strcmp(name,"Pressure")==0) return PressureEnum;
    507               else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
     511              if (strcmp(name,"Pressure")==0) return PressureEnum;
     512              else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
     513              else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
    512514              else if (strcmp(name,"ResetPenalties")==0) return ResetPenaltiesEnum;
    513515              else if (strcmp(name,"SegmentOnIceShelf")==0) return SegmentOnIceShelfEnum;
    … …  
    627629              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    628630              else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    629               else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    630               else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"MinVx")==0) return MinVxEnum;
     634              if (strcmp(name,"MinVel")==0) return MinVelEnum;
     635              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
     636              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    635637              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    636638              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    … …  
    717719              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    718720              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
     721              else if (strcmp(name,"Seaiceatm")==0) return SeaiceatmEnum;
     722              else if (strcmp(name,"Seaiceocean")==0) return SeaiceoceanEnum;
     723              else if (strcmp(name,"SeaiceThickness")==0) return SeaiceThicknessEnum;
     724              else if (strcmp(name,"SeaiceConcentration")==0) return SeaiceConcentrationEnum;
     725              else if (strcmp(name,"SeaiceSpcvx")==0) return SeaiceSpcvxEnum;
     726              else if (strcmp(name,"SeaiceSpcvy")==0) return SeaiceSpcvyEnum;
     727              else if (strcmp(name,"BasalforcingsRhoOcean")==0) return BasalforcingsRhoOceanEnum;
     728              else if (strcmp(name,"BasalforcingsOceanCoef")==0) return BasalforcingsOceanCoefEnum;
     729              else if (strcmp(name,"BasalforcingsOceanLinDragCoef")==0) return BasalforcingsOceanLinDragCoefEnum;
     730              else if (strcmp(name,"BasalforcingsOceanQuadDragCoef")==0) return BasalforcingsOceanQuadDragCoefEnum;
     731              else if (strcmp(name,"BasalforcingsOceanTurningAngle")==0) return BasalforcingsOceanTurningAngleEnum;
     732              else if (strcmp(name,"BasalforcingsOceanSsh")==0) return BasalforcingsOceanSshEnum;
     733              else if (strcmp(name,"BasalforcingsOceanVx")==0) return BasalforcingsOceanVxEnum;
     734              else if (strcmp(name,"BasalforcingsOceanVy")==0) return BasalforcingsOceanVyEnum;
     735              else if (strcmp(name,"SurfaceforcingsRhoAir")==0) return SurfaceforcingsRhoAirEnum;
     736              else if (strcmp(name,"SurfaceforcingsAirCoef")==0) return SurfaceforcingsAirCoefEnum;
     737              else if (strcmp(name,"SurfaceforcingsAirLinDragCoef")==0) return SurfaceforcingsAirLinDragCoefEnum;
     738              else if (strcmp(name,"SurfaceforcingsAirQuadDragCoef")==0) return SurfaceforcingsAirQuadDragCoefEnum;
     739              else if (strcmp(name,"SurfaceforcingsWindVx")==0) return SurfaceforcingsWindVxEnum;
     740              else if (strcmp(name,"SurfaceforcingsWindVy")==0) return SurfaceforcingsWindVyEnum;
     741              else if (strcmp(name,"Matseaice")==0) return MatseaiceEnum;
     742              else if (strcmp(name,"MaterialsPoisson")==0) return MaterialsPoissonEnum;
     743              else if (strcmp(name,"MaterialsYoungModulus")==0) return MaterialsYoungModulusEnum;
     744              else if (strcmp(name,"MaterialsDamage")==0) return MaterialsDamageEnum;
     745              else if (strcmp(name,"MaterialsRidgingExponent")==0) return MaterialsRidgingExponentEnum;
     746              else if (strcmp(name,"MaterialsCohesion")==0) return MaterialsCohesionEnum;
     747              else if (strcmp(name,"MaterialsInternalFrictionCoef")==0) return MaterialsInternalFrictionCoefEnum;
     748              else if (strcmp(name,"MaterialsCompressionCoef")==0) return MaterialsCompressionCoefEnum;
     749              else if (strcmp(name,"MaterialsTractionCoef")==0) return MaterialsTractionCoefEnum;
     750              else if (strcmp(name,"VxStar")==0) return VxStarEnum;
     751              else if (strcmp(name,"VyStar")==0) return VyStarEnum;
    719752              else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
    720753         else stage=7;

  • issm/trunk-jpl/src/m/classes/geometry.m

    r17931 r18492  
    5858                function md = checkconsistency(obj,md,solution,analyses) % {{{
    5959
     60                        if (solution==SeaiceSolutionEnum()),
     61                                return;
     62                        end
     63
    6064                        if (solution==TransientSolutionEnum() & md.transient.isgia) | (solution==GiaSolutionEnum()),
    6165                                md = checkfield(md,'fieldname','geometry.thickness','forcing',1,'NaN',1,'>=',0);

  • issm/trunk-jpl/src/m/classes/mask.m

    r18068 r18492  
    4040                end % }}}
    4141                function md = checkconsistency(obj,md,solution,analyses) % {{{
     42                        if (solution==SeaiceSolutionEnum()),
     43                                return;
     44                        end
    4245
    4346                        md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices 1]);
    … …  
    6164                        WriteData(fid,'object',obj,'fieldname','ice_levelset','format','DoubleMat','mattype',1);
    6265
     66                        if(md.private.solution==SeaiceSolutionEnum()),
     67                                return;
     68                        end
     69
    6370                        % get mask of vertices of elements with ice
    6471                        isice=md.mask.ice_levelset<=0.;

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

    r18480 r18492  
    3939                gia              = 0;
    4040
     41                seaice           = 0;
     42
    4143                autodiff         = 0;
    4244                flaim            = 0;
    … …  
    132134                                case 0
    133135                                        md=setdefaultparameters(md);
     136                                case 1
     137                                        if strcmpi(varargin{1},'seaice'),
     138                                                md=setdefaultparameters(md);
     139                                                md.materials = matseaice();
     140                                                md.surfaceforcings = seaiceatm();
     141                                                md.basalforcings   = seaiceocean();
     142                                                md.initialization  = seaiceinitialization();
     143                                        else
     144                                                error('model constructor not supported yet');
     145                                        end
     146
    134147                                otherwise
    135148                                        error('model constructor error message: 0 of 1 argument only in input.');
    … …  
    11181131                        md.transient        = transient();
    11191132                        md.gia              = gia();
     1133                        md.seaice           = seaice();
    11201134                        md.autodiff         = autodiff();
    11211135                        md.flaim            = flaim();
    … …  
    12951309                        disp(sprintf('%19s: %-22s -- %s','radaroverlay'    ,['[1x1 ' class(obj.radaroverlay) ']'],'radar image for plot overlay'));
    12961310                        disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(obj.miscellaneous) ']'],'miscellaneous fields'));
     1311                        disp(sprintf('%19s: %-22s -- %s','seaice'       ,['[1x1 ' class(obj.seaice) ']'],'parameters for Sea Ice solution'));
    12971312                end % }}}
    12981313                function memory(obj) % {{{

  • issm/trunk-jpl/src/m/classes/stressbalance.m

    r17931 r18492  
    210210                function marshall(obj,md,fid) % {{{
    211211
     212                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
     213
     214                        if md.private.solution==SeaiceSolutionEnum,
     215                                return;
     216                        end
     217
    212218                        yts=365.0*24.0*3600.0;
    213219
    … …  
    223229                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','maxiter','format','Integer');
    224230                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','shelf_dampening','format','Integer');
    225                         WriteData(fid,'object',obj,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
    226231                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','penalty_factor','format','Double');
    227232                        WriteData(fid,'object',obj,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer');

  • issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.m

    r17831 r18492  
    9191                analyses=[DamageEvolutionAnalysisEnum()];
    9292
     93        case SeaiceSolutionEnum(),
     94                analyses=[SeaiceAnalysisEnum()];
     95
    9396        otherwise
    9497                error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!');

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r18476 r18492  
    382382def GiaSolutionEnum(): return StringToEnum("GiaSolution")[0]
    383383def GiaAnalysisEnum(): return StringToEnum("GiaAnalysis")[0]
     384def SeaiceSolutionEnum(): return StringToEnum("SeaiceSolution")[0]
     385def SeaiceAnalysisEnum(): return StringToEnum("SeaiceAnalysis")[0]
    384386def MeshdeformationSolutionEnum(): return StringToEnum("MeshdeformationSolution")[0]
    385387def MeshdeformationAnalysisEnum(): return StringToEnum("MeshdeformationAnalysis")[0]
    … …  
    694696def LevelsetfunctionSlopeYEnum(): return StringToEnum("LevelsetfunctionSlopeY")[0]
    695697def LevelsetfunctionPicardEnum(): return StringToEnum("LevelsetfunctionPicard")[0]
     698def SeaiceatmEnum(): return StringToEnum("Seaiceatm")[0]
     699def SeaiceoceanEnum(): return StringToEnum("Seaiceocean")[0]
     700def SeaiceThicknessEnum(): return StringToEnum("SeaiceThickness")[0]
     701def SeaiceConcentrationEnum(): return StringToEnum("SeaiceConcentration")[0]
     702def SeaiceSpcvxEnum(): return StringToEnum("SeaiceSpcvx")[0]
     703def SeaiceSpcvyEnum(): return StringToEnum("SeaiceSpcvy")[0]
     704def BasalforcingsRhoOceanEnum(): return StringToEnum("BasalforcingsRhoOcean")[0]
     705def BasalforcingsOceanCoefEnum(): return StringToEnum("BasalforcingsOceanCoef")[0]
     706def BasalforcingsOceanLinDragCoefEnum(): return StringToEnum("BasalforcingsOceanLinDragCoef")[0]
     707def BasalforcingsOceanQuadDragCoefEnum(): return StringToEnum("BasalforcingsOceanQuadDragCoef")[0]
     708def BasalforcingsOceanTurningAngleEnum(): return StringToEnum("BasalforcingsOceanTurningAngle")[0]
     709def BasalforcingsOceanSshEnum(): return StringToEnum("BasalforcingsOceanSsh")[0]
     710def BasalforcingsOceanVxEnum(): return StringToEnum("BasalforcingsOceanVx")[0]
     711def BasalforcingsOceanVyEnum(): return StringToEnum("BasalforcingsOceanVy")[0]
     712def SurfaceforcingsRhoAirEnum(): return StringToEnum("SurfaceforcingsRhoAir")[0]
     713def SurfaceforcingsAirCoefEnum(): return StringToEnum("SurfaceforcingsAirCoef")[0]
     714def SurfaceforcingsAirLinDragCoefEnum(): return StringToEnum("SurfaceforcingsAirLinDragCoef")[0]
     715def SurfaceforcingsAirQuadDragCoefEnum(): return StringToEnum("SurfaceforcingsAirQuadDragCoef")[0]
     716def SurfaceforcingsWindVxEnum(): return StringToEnum("SurfaceforcingsWindVx")[0]
     717def SurfaceforcingsWindVyEnum(): return StringToEnum("SurfaceforcingsWindVy")[0]
     718def MatseaiceEnum(): return StringToEnum("Matseaice")[0]
     719def MaterialsPoissonEnum(): return StringToEnum("MaterialsPoisson")[0]
     720def MaterialsYoungModulusEnum(): return StringToEnum("MaterialsYoungModulus")[0]
     721def MaterialsDamageEnum(): return StringToEnum("MaterialsDamage")[0]
     722def MaterialsRidgingExponentEnum(): return StringToEnum("MaterialsRidgingExponent")[0]
     723def MaterialsCohesionEnum(): return StringToEnum("MaterialsCohesion")[0]
     724def MaterialsInternalFrictionCoefEnum(): return StringToEnum("MaterialsInternalFrictionCoef")[0]
     725def MaterialsCompressionCoefEnum(): return StringToEnum("MaterialsCompressionCoef")[0]
     726def MaterialsTractionCoefEnum(): return StringToEnum("MaterialsTractionCoef")[0]
     727def VxStarEnum(): return StringToEnum("VxStar")[0]
     728def VyStarEnum(): return StringToEnum("VyStar")[0]
    696729def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]
Note: See TracChangeset for help on using the changeset viewer.