Wiki
Clone wikigeomIO / VaryGeomTutorial
[Home] [Basic 2D tutorial] [Basic 3D tutorial] [VaryGeom2D tutorial] [VaryGeom3D tutorial]
Example
As with all geomIO functions, make sure that the geomIO directory is added to your Matlab path. You can reproduce the following example by using the files located in geomio/tutorials/TutoVaryGeom2D/.
Preparation
We will work with a simple 2D example of a salt diapir given in Dipir.EW.svg
as it makes many of the following illustrations more intuitive. We start with some basic geomIO options:
% set rng seed for reproducibility rng(200877) % basic geomIO settings inputFile = './Input/Diapir.EW.svg'; opt = geomIO_Options(); opt.outputDir = './Output_Salt'; opt.inputFileName = inputFile; opt.LaMEMinputFileName = './Input/Salt.dat'; opt.readLaMEM = true; opt.writeParaview = true; opt.writePolygons = true; opt.interp = true; opt.zi = [-9:1:9]; opt.getPolygons = true;
In the next step, the different volumes present in the drawing are defined:
% volumes paths = { 'Air', 0, 0 'Crust', 1, 0 'Diapir', 2, 0 'SaltBed', 2, 0 }; opt.pathNames = {paths{:,1}}; opt.phase = [paths{:,2}]; opt.type = [paths{:,3}];
From here, geomIO can be run to generate the reference model and visualize it to pick appropriate control polygons and scaling parameters. The following code snippet also produces .vtk files of the volumes to view in Paraview (opt.writeParaview = true
) and a polygon file containing the entire setup to be read by LaMEM (opt.writePolygons = true
).
% run geomIO [PathCoord,Volumes,opt] = run_geomIO(opt,'default'); % visualize diapir polygons = Volumes.Diapir.Polygons; figure() hold on; box on for i = 1 : length(polygons) plot(polygons{i}(:,1),polygons{i}(:,3),'k') end
% visualize control polygons CtrlPoly = [153, 253, 385, 482]; for i = 1 : 4 plot(polygons{CtrlPoly(i)}(:,1),polygons{CtrlPoly(i)}(:,3),'r','LineWidth',2) end
Vary geometry
With the control polygons picked, we can set up the options for variable geometry. For a description of the individual options, check the list given at the bottom of the page.
% vary geom settings opt.varyGeom.do = true; % activate geometry vatriation opt.varyGeom.volume = 'Diapir'; % name of volume in Inkscape opt.varyGeom.dim = '2D_X'; % mode opt.varyGeom.genType = 'grid'; % how to generate scaling parameters opt.varyGeom.Grid.lb = [0.5; 0.5; 0.5; 0.5]; % lower bound for scaling parameter generation opt.varyGeom.Grid.ub = [2.0; 3.0; 2.0; 1.5]; % upper bound for scaling parameter generation opt.varyGeom.Grid.num = [2; 2; 2; 2]; % number of samples for each scaling parameter opt.varyGeom.Grid.Sz = [0.8; 1.1; 2]; % options for vertical scaling parameter generation opt.varyGeom.Grid.noise = 0.3; % noise to be added to scaling parameters opt.varyGeom.stretchType = 'factor'; % scale by factor or absolute value opt.varyGeom.SzType = 'bot'; % reference level for vertical stretching opt.varyGeom.CtrlPoly = [153, 253, 385, 482]; % indices of control polygons opt.varyGeom.CP_ref = 'global'; % reference system for control polygons opt.varyGeom.CP_track = true; % control polygons cahnge depth with vertical scaling opt.varyGeom.outName = './SaltExample'; % output name opt.varyGeom.outNameOffset = 0; % offset numbering of output opt.varyGeom.writeParaview = false; % write .vtk files opt.varyGeom.writePolygons = true; % write LaMEM polygons opt.varyGeom.drawOutlines = true; % draw outline along profiles opt.varyGeom.outlineProf = [0, 0]; % coordinates of outline profiles
Running geomIO with these additional options produces 32 different setups for LaMEM to be read (opt.varyGeom.writePolygons = true
) with identical background lithosphere and different salt diapirs. The scaling parameters are generated on a regular grid within the boundaries provided by the options opt.varyGeom.Grid.lb
and opt.varyGeom.Grid.ub
.
% run geomIO [PathCoord,Volumes,opt] = run_geomIO(opt,'default');
opt.varyGeom.drawOutlines = true
, geomIO also outputs a .mat file containing profiles of the variations it produced. They can be loaded and visualized with the following code snippet, showing the reference geometry in red.
% visualize produced variations load Output_Salt/Outlines.mat figure() hold on, box on for i = 1 : size(Outline,2)-1 plot(Outline{1,i}(:,1),Outline{1,i}(:,2),'k') end plot(Outline{1,33}(:,1),Outline{1,33}(:,2),'r','LineWidth',2) axis([-10 10 -16.5 0])
Note that in the .svg file, we also included a salt layer extending from 16.5 km depth downwards, overprinting the deepest part of the diapir for all variations. There are gaps between the variations because the scaling parameters were generated on a regular grid. By switching to a different generation mode, the scaling parameters can instead be generated by using normal distributions, resulting in a more continuous result.
% opt.varyGeom.genType = 'grid'; % how to generate scaling parameters % opt.varyGeom.Grid.lb = [0.5; 0.5; 0.5; 0.5]; % lower bound for scaling parameter generation % opt.varyGeom.Grid.ub = [2.0; 3.0; 2.0; 1.5]; % upper bound for scaling parameter generation % opt.varyGeom.Grid.num = [2; 2; 2; 2]; % number of samples for each scaling parameter % opt.varyGeom.Grid.Sz = [0.8; 1.1; 2]; % options for vertical scaling parameter generation % opt.varyGeom.Grid.noise = 0.3; % noise to be added to scaling parameters opt.varyGeom.genType = 'normal'; % how to generate scaling parameters opt.varyGeom.n = 32; % number of variations opt.varyGeom.mu = 1; % mean of scaling parameters opt.varyGeom.sigma = 0.3; % standard deviation of scaling parameters opt.varyGeom.Sz = [1, 0.08]; % mean and standard deviation of vertical scaling
geomIO will ouput the scaling parameters it generated as Output_Salt/StretchInfo.mat
. You can also generate the scaling parameters yourself and instead of using opt.varyGeom.genType
, utilize opt.varyGeom.loadStretch = yourScalingParams.mat
. Make sure the file matches the format of Output_Salt/StretchInfo.mat
.
Subduction zones
For a subduction zone setting, the basic geomIO settings and way, volumes are read from the .svg file are identical. Here, we use the Subduction.EW.svg
as input.
% basic geomIO settings inputFile = './Input/Subduction.EW.svg'; opt = geomIO_Options(); opt.outputDir = './Output_Subduction'; opt.inputFileName = inputFile; opt.LaMEMinputFileName = './Input/Subduction.dat'; opt.readLaMEM = true; opt.writeParaview = true; opt.writePolygons = true; opt.interp = true; opt.zi = [-9:1:9]; opt.getPolygons = true; % volumes paths = { 'Mantle', 6, 0 'ContCrust', 1, 0 'ContMantle', 2, 0 'SimpleSlab', 3, 0 }; opt.pathNames = {paths{:,1}}; opt.phase = [paths{:,2}]; opt.type = [paths{:,3}];
To change the angle of the slab, we use a different set of options. With these options, we can place rotation centers along the subducting plate. At these locations, the part of the plate opposite the trench is rotated to steepen or flatten it. The following options steepen the slab at the trench by 10° and and flatten it by 20° 120 km away from the trench.
% vary subduction settings opt.varySub.do = true; % activate subduction variation opt.varySub.vols = 'SimpleSlab'; % name of volume in Inkscape opt.varySub.ref = 'trench'; % reference system of rotation center coordinates opt.varySub.xRot = [0, 120]; % rotation center coordinates opt.varySub.theta = [10, -20]; % angles of rotation opt.varySub.tolZ = 1; % tolerance when identifying the plate opt.varySub.drawOutlines = true; % draw outline along profiles opt.varySub.outlineProf = [0, 0]; % coordinates of outline profiles
% options to add layers opt.varySub.addWZ = true; % add a weak zone opt.varySub.d_WZ = 20; % thickness of weak zone opt.varySub.ID_WZ = 4; % LaMEM ID opt.varySub.type_WZ = 0; % LaMEM marker type opt.varySub.d_Lith = 100; % thickness of overriding plate opt.varySub.addCrust = true; % add crust opt.varySub.d_OC = 10; % thickness of crust opt.varySub.ID_OC = 5; % LaMEM ID opt.varySub.type_OC = 0; % LaMEM marker type
d_WZ
and d_OC
control the thickness of the layers respectively, and d_Lith controls how deep the weak zone goes. This should fit the thickness of the overrding plate. Finally, just launch geomIO as always.
% Run geomIO [PathCoord,Volumes,opt] = run_geomIO(opt,'default');
To visualize the result, use the following code.
% Visualize load ./Input/OGSub_Outline.mat load Output_Subduction/Sub_Outline.mat figure() hold on; box on plot(OG_Outlines.Plate.EW(:,1),OG_Outlines.Plate.EW(:,2),'k:','LineWidth',0.5) plot(SubOutlines.Plate.EW(:,1),SubOutlines.Plate.EW(:,2),'k') plot(SubOutlines.OC.EW(:,1),SubOutlines.OC.EW(:,2),'b') plot(SubOutlines.WZ.EW(:,1),SubOutlines.WZ.EW(:,2),'r') axis equal axis([-350 200 -350 50])
Dotted black line shows the original plate drawn in Inkscape, solid black line shows the new variation, blue line shows the oceanic crust that was automaticall generated and red line the weak zone that was automatically generated.
All options
In the following, all available options are listed. They are split into different categories here but the order of use does not matter. To use an option, add it to the matlab script as opt.varyGeom.do = true
. You don't need to set all options. Flag options are defaulted to false
. If an important option is missing, the code will throw errors poiting out the missing options.
General options
Option | Description | Accepted input |
---|---|---|
do | Flag that activates the feature. | true, false |
volume | Name of volume to vary. | a string |
n | Number of variations to generate | an integer |
stretchType | Scale by factor (default) or absolute value. | 'factor', 'absolute' |
dim | Which dimension to vary. 3D: uses Sx and Sy. 3D_homo: uses Sx=Sy. 2D_X: uses Sx and Sy=1. 2D_Y uses Sy and Sx=1. | '3D', '3D_homo', '2D_X', '2D_Y' |
Control polgons and scaling parameters
Option | Description | Accepted input |
---|---|---|
CtrlPoly | Indices of polygons that are control polygons. | ordered row vector of integers (e.g. [3,7,17]) |
CP_ref | Indices will be in reference to what. global(default): entire model. volume: the volume in question. | 'global', 'volume' |
CP_track | Flag that activates shfiting control polygons when volume is scaled in z-direction. | true, false |
SzType | Which z-coordinate of the volume stays fixed? com: center of mass. top: top of the volume. bottom: bottom of the volume. middle: center of the volume. | 'com', 'top', 'bottom', 'middle' |
Smin | Minimum scaling factor to prevent polygons becoming too small. Default is 0.01. | a float |
loadStretch | Load your own set of scaling parameters instead of making geomIO generate them, name of .mat file to load. | a string |
genType | Have geomIO generate the scaling parameters. normal: uses a normal distribution. grid: generates parameters on a regular grid. chain: uses a markov-chain. guidedChain: uses a guided markow-chain. | 'normal', 'grid', 'chain', 'guidedChain' |
mu | Required for genType: normal, chain, guidedChain. Mean of Sx and Sy. | one or two floats depedning on dim |
sigma | Required for genType: normal, chain, guidedChain. Standard deviation of Sx and Sy. | one or two floats depending on dim |
sigma2 | Option for genType: chain, guidedChain. Controls how close the next values is to the previous one. Default is sigma. | a float |
Sz | Option for genType: normal, chain, guidedChain. Use if scaling in z-direction is desired. Mean of Sz and sigma of Sz. | two floats |
Grid.lb | Required for genType: grid. Lower bound of values. | column vector, one row for each control polygon, 2 columns for dim='3D' |
Grid.ub | Required for genType: grid. Upper bound of values. | column vector, one row for each control polygon, 2 columns for dim='3D' |
Grid.num | Required for genType: grid. Number of values to be selected between lb and ub. | column vector, one row for each control polygon, 2 columns for dim='3D' |
Grid.Sz | Option for genType: grid. lb, ub and num for Sz if scaling in z-direction is desired. | [float, float, integer] |
Grid.noise | Option for genType: grid. Adds random noise to parameters. Value will be multiplied with half of grid spacing and serve as standard deviation. Should be < 0.5. | a float |
Output options
Option | Description | Accepted input |
---|---|---|
outName | Name of output files. | a string |
outNameOffset | A continous number is added to output files. To add to output of previous uses, the numbering can start at an desired value. | an integer |
writeParaview | Flag that activates writing a .vtk files for each variation (a few MB for each variation and volume). | true, false |
writePolygons | Flag that activates writing a polygon file to be read as inout by LaMEM for each variation (a few MB for each variation). | true, false |
saveVols | Save the produced volumes. 0: No. 1: Yes, but only strictly necessary information. 2: Yes, save everything. | 0, 1, 2 |
drawOutlines | Flag that outputs a file containing the outlines of each variation along a EW and NS profile. | true, false |
outlineProf | Y-coordinate of EW-profile and X-ccordinate of NS-profile. | [float, float] |
Additional options
Option | Description | Accepted input |
---|---|---|
loadPolygons | Replace a volume from the .svg with one that is loaded from a .mat file. | ['path/to/mat/file.mat', 'nameOfVolume'] |
gravOnly | Flag to activate a simpler algorithm for scaling in z-direction. Not LaMEM compatible. Can be used with geomIO's gravity forward modeling. | true, false |
shift | Adds a shift of the entire body in the x-y-plane. Mean of shift_x and shift_y and standard deviation of shift_x and shift_y. | [float, float, float, float] |
shiftZ | Adds a shift of the entire body in z-direction. Mean and sigma of shift_z | [float, float] |
theta | Rotates coordinate system clockwise before scaling is applied. | a float |
All subduction options
Subduction zones are handled differently internally and have their own set of options listed below. To use an option, add it to the matlab script as opt.varySub.do = true
. All flag options are defaulted to 'false'.
Basic options
Option | Description | Accepted input |
---|---|---|
do | Flag that activates the feature. | true, false |
vols | Name of subducting plate volume. Currently only one plate is possible. | a string |
ref | Reference of the coordinate system for following options. coord (default): uses absolute coordinates. trench: uses coordinates relative to trench location. | 'coord', 'trench' |
xRot | Distance of rotation center from reference. Provide one entry per rotation center. This will be parallel to drawing profiles. | ordered row vectorof floats (e.g., [0, 150.5, 300]) |
theta | Change in dip angle for everything beyond xRot (>0: steepen, <0: flatten). Provide one entry per rotation center. | ordered row vector of floats (e.g., [15, -12.5, -10]) |
tolZ | Tolerance in z-direction when detecting the flat part of the slab. Default is 2. | a float |
drawOutlines | Flag that outputs a file containing the outlines of plate, crust and weak zone along a EW and NS profile. | true, false |
outlineProf | Y-coordinate of EW-profile and X-ccordinate of NS-profile. | [float, float] |
Add layers
Option | Description | Accepted input |
---|---|---|
addWZ | Flag to add a weak zone on top of the slab. | true, false |
d_WZ | Thickness of the weak zone. | a float |
ID_WZ | Phase ID of the weak zone to match in LaMEM. | an integer |
type_WZ | Marker type of weak zone in LaMEM. 0 (default):overwriting. 1: additive. | 0, 1 |
d_Lith | Depth to which the weak zone should extend. This should match the thickness of the overriding plate. | a float |
addCrust | Flag to turn the top layer of the plate into crust. So you don't have to draw it individually. | true, false |
d_OC | Thickness of the crust. | a float |
ID_OC | Phase ID of the crust to match in LaMEM. | an integer |
type_OC | Marker type of the crust in LaMEM. 0 (default):overwriting. 1: additive. | 0, 1 |
Questions/Suggestions
Send questions and suggestions about the tool to arspang@uni-mainz.de.
How to cite
If you use any functionality described here, please cite:
Bauville, A., & Baumann, T. S.(2019). geomIO: An open-source MATLAB toolbox to create the initial configuration of 2-D/3-D thermo-mechanical simulations from 2-D vector drawings. Geochemistry, Geophysics, Geosystems, 20, 1665–1675. https://doi.org/10.1029/2018GC008057
and
Spang, A., Baumann, T. S., & Kaus, B. J. P. (2022). Geodynamic modeling with uncertain initial geometries. Geochemistry, Geophysics, Geosystems, 23, e2021GC010265. https://doi.org/10.1029/2021GC010265
Updated