[Home] [Basic 2D tutorial] [Basic 3D tutorial] [VaryGeom2D tutorial] [VaryGeom3D tutorial]

geomIO_logo

Example
- Preparation
- Vary geometry
Subduction zones
All options
All subduction options
- Basic options
- Add layers
Questions/Suggestions
How to cite

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

In this example, we pick 4 control polygons (153, 253, 385, 482) that represent the diapir's base, neck, transition from neck to head and head. They are marked red.

% 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');

Because of the option 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

A weak zone on top of the slab and an oceanic crust within the slab can easily added with the following options.

% 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

Options 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

Wiki

geomIO / VaryGeomTutorial