[26500] | 1 | # setmask {{{
|
---|
| 2 | """
|
---|
| 3 | SETMASK - establish boundaries between grounded and floating ice.
|
---|
| 4 |
|
---|
| 5 | By default, ice is considered grounded. The contour floatingicename defines nodes
|
---|
| 6 | for which ice is floating. The contour groundedicename defines nodes inside an floatingice,
|
---|
| 7 | that are grounded (ie: ice rises, islands, etc ...)
|
---|
| 8 | All input files are in the Argus format (extension .exp).
|
---|
| 9 |
|
---|
| 10 | Usage:
|
---|
| 11 | md=setmask(md,floatingicename,groundedicename)
|
---|
| 12 |
|
---|
| 13 | Examples:
|
---|
| 14 | md=setmask(md,'all','');
|
---|
| 15 | md=setmask(md,'Iceshelves.exp','Islands.exp');
|
---|
| 16 | """
|
---|
| 17 | function setmask(md::model,floatingicename::String,groundedicename::String)
|
---|
| 18 |
|
---|
[26713] | 19 | elementonfloatingice = FlagElements( md, floatingicename)
|
---|
| 20 | elementongroundedice = FlagElements( md, groundedicename)
|
---|
[26500] | 21 |
|
---|
| 22 | elementonfloatingice = convert( Array{Float64}, (elementonfloatingice.>0) .& (elementongroundedice.==0.))
|
---|
| 23 | elementongroundedice = convert( Array{Float64}, elementonfloatingice.==0.)
|
---|
| 24 |
|
---|
| 25 | vertexonfloatingice=zeros(md.mesh.numberofvertices)
|
---|
| 26 | vertexongroundedice=zeros(md.mesh.numberofvertices)
|
---|
| 27 |
|
---|
| 28 | vertexongroundedice[md.mesh.elements[findall(elementongroundedice.>0),:]] .= 1.
|
---|
| 29 | vertexonfloatingice[findall(vertexongroundedice.==0.)] .= 1.
|
---|
| 30 |
|
---|
| 31 | #define levelsets
|
---|
| 32 | md.mask.ocean_levelset = vertexongroundedice
|
---|
| 33 | md.mask.ocean_levelset[findall(vertexongroundedice .==0.)] .= -1.
|
---|
| 34 | md.mask.ice_levelset = -1*ones(md.mesh.numberofvertices)
|
---|
| 35 |
|
---|
| 36 | return md
|
---|
| 37 | end
|
---|
| 38 | #}}}
|
---|
| 39 | # FlagElements {{{
|
---|
| 40 | function FlagElements(md::model,region::String)
|
---|
| 41 |
|
---|
| 42 | if isempty(region)
|
---|
| 43 | flags = zeros(md.mesh.numberofelements)
|
---|
| 44 | elseif region == "all"
|
---|
| 45 | flags = ones(md.mesh.numberofelements)
|
---|
| 46 | else
|
---|
[26713] | 47 | xcenter = md.mesh.x[md.mesh.elements]*[1;1;1]/3
|
---|
| 48 | ycenter = md.mesh.y[md.mesh.elements]*[1;1;1]/3
|
---|
| 49 | flags = ContourToNodes(xcenter, ycenter, region, 2.)
|
---|
[26500] | 50 | end
|
---|
| 51 |
|
---|
| 52 | return flags
|
---|
| 53 | end
|
---|
| 54 | #}}}
|
---|