-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfric2d_workflow.m
1498 lines (1129 loc) · 55.7 KB
/
fric2d_workflow.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
function GROWOutputStruct = fric2d_workflow(fileName,X,Y,varargin)
% Automated General work flow for the fric2d software
% INPUT:
% filename: string specifying the name of the fric2D input file (no
% extension)
% X: N by 1 array of x-coordinates of a desired fault profile
% Y: N by 1 array of y-coordinates of a desired fault profile
%
% varargin:
% not yet set...
% OUTPUT:
%
% GROWOUTPUT = fric2d_workflow(fileName,X,Y) - (if one output)
% GROWOUTPUT: structure with with fields corresponding to blocks of data
% segments: 1 by N (N is numuber of fault segments, or, equivalently, the
% number of defects+1) structure with fields corresponding to .eff
% headers for the crack information (num xbeg ybeg xend yend
% stiffS stiffN ten_str init_coh slid_coh stat_fric dy_fric crit_slip)
%
% must be in same folder as:
% - simple_shear_bc.pl
% - fric2D files
% this file create the boundary line geometry and boundary stress
% conditions (this might be changed in future edits so that all inputs
% are set within the matlab environment - specifically for rock properties)
% must have the following paths to programs
% to fric2d (version fric2d_3.2.8): fric2d_code/fric2d_3.2.8 -i
% to GROW (version GROW_JULY_2016.pl): ./GROW_JULY_2016.pl
%--------------------------------------------------------------------------
% system commands are for set for usage with mac or linux
% -------------------------------------------------------------------------
%% input parameters:
% -----------------------------------------------------------------------
[inputParameters, GROWInputParameters] = userInput();
% -----------------------------------------------------------------------
%% work flow:
% *check Input*:
[X, Y,plotOption] = checkinput(fileName, X,Y,varargin, nargin);
% *check Directory*
check_directory(GROWInputParameters);
%% 1) make input file
[~, newX] = makeinputfile(fileName,X,Y,inputParameters);
%% 2) run program
runfric2D(fileName);
%% 3) extract info
outputStruct = extractoutput(fileName);
%% 4) assess shear failure
[shearFailureCoord , tensileFailureCoord, ...
coulombCriterion , tensileCriterion] = assessfailure(newX,Y, outputStruct, inputParameters);
%% 5) feed into GROW
% 5a)choose which crack(s) to grow (now set to 'all', could alternatively
% be one of: 'all', 'tensile', 'shear', 'max coulomb', and 'max tensile')
crackLocations = choosecracks(shearFailureCoord , tensileFailureCoord , ...
inputParameters.user.cracks2grow, inputParameters);
if ~isempty(crackLocations)
% 5b) run GROW with the chosen ckracks
%runGROW(fileName, crackLocations, meanSegmentLength, inputParameters, GROWInputParameters);
runGROW(newX, Y, crackLocations,inputParameters, GROWInputParameters);
%% 6) extract GROW output (not finished)
GROWOutputStruct = extractGROWoutput(crackLocations, GROWInputParameters.growFileName, inputParameters);
end
%% 7) graphically show entire output
plotoutputinfo;
%% 8) organize files
% create directory where all files will be stored
organizefiles;
%% 9) check model - issue potential warnings:
checkmodel(X,Y,inputParameters, outputStruct, GROWOutputStruct);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ----------------------- embeded functions ------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 7) plot output information into three plots
function [] = plotoutputinfo(varargin)
if ~strcmp(inputParameters.user.showFric2DOutput,'yes')
if strcmp(inputParameters.user.showFric2DOutput, 'no')
disp('output not displayed change inputParameters.user.showFric2DOutput in input section to ''yes'' if so desired')
else
error('inputParameters.user.showFric2DOutput in input section must be ''yes'' or ''no''')
end
return
end
% three plots:
% a) geometry
% b) displacements
% c) stresses
% d) failure criterion
%%%%%%%%%%%%%%%%%%%%%%%%%
figure
longPlotOn = strcmp(plotOption,'long plot');
numPlots = 4;
subPlotCount = 1;
if longPlotOn
subplot(numPlots, 1, subPlotCount)
end
subPlotCount = subPlotCount+1;
% a) geometry: (boundary box, fault, failure points)
hold on
% plot the poundary box as a rectangle
boundaryBox = fill( outputStruct.boundaryLines.xBeg , ...
outputStruct.boundaryLines.yBeg , ...
[0.8 0.8 0.8]);
% make quiver plot of boundary conditions
angle = outputStruct.boundaryLines.angle;
US = outputStruct.boundaryLines.USMonotonic;
UN = outputStruct.boundaryLines.USMonotonic;
hold on
% make vectors
elementDir =[cosd(angle),sind(angle)]./repmat((sqrt(cosd(angle).^2+sind(angle).^2)),1,2); %unit vecto in the direction of the element;
USVec = elementDir.*repmat(US,1,2);
R = [0 -1; 1, 0];
UNVec = (R*elementDir'.*repmat(UN,1,2)')'; % rotated 90 degrees
USQuivers = quiver(outputStruct.boundaryLines.xBeg , outputStruct.boundaryLines.yBeg , ...
USVec(:,1), USVec(:,2), ...
'AutoScale', 'off');
UNQuivers1 = quiver(outputStruct.boundaryLines.xBeg , outputStruct.boundaryLines.yBeg , ...
UNVec(:,1), UNVec(:,2), ...
'AutoScale', 'off');
UNQuivers2 = quiver(outputStruct.boundaryLines.xBeg , outputStruct.boundaryLines.yBeg , ...
-UNVec(:,1), -UNVec(:,2), ...
'AutoScale', 'off');
% plot the fault segments outputed from grow (if clause is within the function)
plotGROWOutput('geometry')
% superpose the fault geometry
faultProfile = plot( outputStruct.faultInfo.xBeg , ...
outputStruct.faultInfo.yBeg , ...
'.-');
% plot endbits
leftEndBitLine = plot( outputStruct.leftEndBit.xBeg, ...
outputStruct.leftEndBit.yBeg, 'r.-');
rightEndBitLine = plot( outputStruct.rightEndBit.xBeg, ...
outputStruct.rightEndBit.yBeg, 'r.-');
% plot the points that a subject to failure (local maxima)
tensileFailurePoints = scatter(tensileFailureCoord(:,1), tensileFailureCoord(:,2));
shearFailurePoints = scatter(shearFailureCoord(:,1), shearFailureCoord(:,2));
hLegendA = legend([boundaryBox , ...
faultProfile , ...
rightEndBitLine , ...
tensileFailurePoints , ...
shearFailurePoints ], ...
...
'Boundary Elements' , ...
'Rough Fault Profile' , ...
'end bits' , ...
'Failure in tension' , ...
'Failure in shear' );
set( gca , ...
'Box' , 'off' , ...
'FontName' , 'Helvetica' );
axis equal
xlabel('x (m)')
ylabel('y (m)')
if longPlotOn
% b) displacmement (shear, normal)
subplot(numPlots, 1, subPlotCount)
subPlotCount = subPlotCount+1;
hold on
yyaxis left
shearDisplacementPlot = plot(outputStruct.faultInfo.xBeg, ...
outputStruct.faultStress.DS);
ylabel('Shear Displacement (m)')
yyaxis right
normalDisplacementPlot = plot(outputStruct.faultInfo.xBeg, ...
outputStruct.faultStress.DN);
ylabel('Normal Displacement (m)')
xlabel('x (m)')
hLegendB = legend([shearDisplacementPlot , ...
normalDisplacementPlot],...
...
'Shear displacement' , ...
'Normal displaceement');
% c) stresses (shear, normal, tangential)
subplot(numPlots, 1, subPlotCount)
subPlotCount = subPlotCount+1;
hold on
shearStressPlot = plot(outputStruct.faultInfo.xBeg, ...
outputStruct.faultStress.SigmaS);
normalStressPlot = plot(outputStruct.faultInfo.xBeg, ...
outputStruct.faultStress.SigmaN);
tangentialStressPlot = plot(outputStruct.faultInfo.xBeg(2:end-1), ...
outputStruct.faultStress2.Tangential);
hLegendC = legend([shearStressPlot , ...
normalStressPlot , ...
tangentialStressPlot],...
...
'Shear Stress' , ...
'Normal Stress' , ...
'Tangential Stress' );
xlabel('x (m)')
ylabel('(MPa)')
% d) failure criterion (coulomb, tensile)
subplot(numPlots, 1, subPlotCount)
subPlotCount = subPlotCount+1;
hold on
coulombPlot = plot(outputStruct.faultInfo.xBeg(2:end-1), ...
coulombCriterion);
tensileStressPlot = plot(outputStruct.faultInfo.xBeg(2:end-1), ...
tensileCriterion);
hLegendD = legend([coulombPlot , ...
tensileStressPlot ], ...
...
'coulomb failure criterion (\tau_m - (\sigma_m sin\phi+c cos\phi) > 0)', ...
'tensile failure criterion (\sigma_3 - c > 0)');
xlabel('x (m)')
ylabel('(MPa)')
end
% nested function:
function plotGROWOutput(what2plot)
% plottting grow output
% input:
% what2plot: cell array, or string with thing to plot
% 'geometry': plot crack geometry
if strcmp(inputParameters.user.runGROWModel, 'yes') % no error meassage here because it should be issued in earlier functions
if ~strcmp(inputParameters.user.plotGROWModel, 'yes')
if strcmp(inputParameters.user.plotGROWModel, 'no')
disp('Grow output not plotted, change inputParameters.user.plotGROWModel to: ''yes''')
else
error('''plotGROWMOdel'' must be ''yes'' or ''no''')
end
return
end
if any(strcmp(what2plot,'geometry'))
hold on
for iSegment = 1:(size(crackLocations,1)+1)
xBeg = GROWOutputStruct.segments(iSegment).xbeg;
yBeg = GROWOutputStruct.segments(iSegment).ybeg;
xEnd = GROWOutputStruct.segments(iSegment).xend;
yEnd = GROWOutputStruct.segments(iSegment).yend;
dx = xEnd - xBeg;
dy = yEnd - yBeg;
segmentLine = quiver(xBeg, yBeg, dx, dy);
segmentLine.ShowArrowHead = 'off';
segmentLine.AutoScale = 'off';
segmentLine.Marker = '.';
end % for segments
else
error('input to function plotGROWOutput must be one, or many, of ''geometry''')
end % if geomtry is called
end % if grow model is run
end % function to plot grow output
end % function to plot output
% 8) organize files to reduce clutter in directory
function organizefiles
% make new directory for the files from this run
% make an output directory if it does not already exist
outputDirectoryName = 'output';
if ~isdir(outputDirectoryName)
system(['mkdir ', outputDirectoryName])
end
newDirectoryName = [outputDirectoryName, '/', fileName];
if ~isdir(newDirectoryName)
command = sprintf('mkdir %s', newDirectoryName);
system(command);
end
% move all file that start with 'fileName' into new directory
command = sprintf('mv %s %s', [fileName,'.*'], [newDirectoryName,'/']);
system(command);
command = sprintf('mv %s %s', 'grow_specs*', [newDirectoryName,'/']);
system(command);
% create input and output directories
if ~isdir([newDirectoryName, '/inputs'])
command = sprintf('mkdir %s', [newDirectoryName, '/inputs']);
system(command);
end
if ~isdir([newDirectoryName, '/outputs'])
command = sprintf('mkdir %s', [newDirectoryName, '/outputs']);
system(command);
end
if ~isdir([newDirectoryName, '/efficient_files'])
command = sprintf('mkdir %s', [newDirectoryName, '/efficient_files']);
system(command);
end
if ~isdir([newDirectoryName, '/prev'])
command = sprintf('mkdir %s', [newDirectoryName, '/prev']);
system(command);
end
% move inpputs an outputs to their respective locations
command = sprintf('mv %s %s', [newDirectoryName,'/*.in'], [newDirectoryName,'/inputs']);
system(command);
command = sprintf('mv %s %s', [newDirectoryName,'/*.out'], [newDirectoryName,'/outputs']);
system(command);
command = sprintf('mv %s %s', [newDirectoryName,'/*.eff'], [newDirectoryName,'/efficient_files']);
system(command);
command = sprintf('mv %s %s', [newDirectoryName,'/*.prev'], [newDirectoryName,'/prev']);
system(command);
% save matlab content including output figure and input specs
savefig([newDirectoryName,'/',fileName]);
save([newDirectoryName,'/',fileName], 'inputParameters', 'GROWInputParameters');
system('rm end_bits.txt');
system('rm input');
system('rm spliced_profile_temp.txt');
system('rm fault_header.txt')
profileName = [fileName,'_profile.txt']; system(['rm ',profileName]);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ------------------- User Input -----------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [inputParameters, GROWInputParameters] = userInput()
inputParameters = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS IS LIKELY THE ONLY SECTION A USER MAY WANT TO CHANGE %%%%%%%%%%%
% material parameters:
E = 3000; % Young's Modulus MPa
nu = 0.25; % Poisson's ratio
T = 10; % Tensile strength of the host rock (MPa) | 5-10
S_0 = 25; % Shear strength of the host rock (MPa) | 25-60 (realistic amounts)
static_friction = 0.6; % static friction of crack elements
dynamic_friction = 0; % dynamic friction of crack elements
critical_slip_distance = 10^-5; % critical slip distance (m)
shear_stiffness = 10^10; % shear stiffness
normal_stiffness = 10^10; % normal stiffness
sliding_cohesion = 0; % sliding cohesion (resistance of fault elements to tension)
% boundary conditions
lithostatic_stress = -40;
shear_stress = -40;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MORE TECHNICAL CHANGES:
% user choices
inputParameters.user.showFric2DOutput = 'yes'; % plot output from fric2d ('yes' or 'no')
inputParameters.user.cracks2grow = 'max coulomb';
% options:
% 'all': all the coordinates specified in input (tensile+shear)
% 'tensile': only cracks that failed in tensile stress (listed in
% tensileFailureCoord)
% 'shear' only cracks that failed in shear stress (listed in
% shearFailureCoord)
% 'max tensile' most tensile crack
% 'max coulomb' most sheared crack
% 'max tensile and max coulomb' both most shear and most tensile
inputParameters.user.runGROWModel = 'yes' ; % run grow ('yes' or 'no')
inputParameters.user.plotGROWModel = 'yes' ; % plot grow output 'yes' or 'no' (make sure runGROWmodel is set to yes as well)
% GROW inputs
GROWInputParameters = [];
GROWInputParameters.growFileName = 'grow_specs';
% grow input parameters parameter
% this will be part of the command to run GROW (Must be digits - no
% decimal points)
GROWInputParameters.angleResolution = 15;
GROWInputParameters.startAngle = 90;
GROWInputParameters.endAngle = 270;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STORING INPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inputParameters.material.E = E; % Young's Modulus MPa
inputParameters.material.nu = nu; % Poisson's ratio
inputParameters.material.c = T; % cohesion (MPa)
inputParameters.material.phi = atand(static_friction); % angle of intenral friction (degrees)
% boundary conditions
inputParameters.boundaryConditions.sigxx = lithostatic_stress; % normal stress parallel to faults MPa
inputParameters.boundaryConditions.sigyy = lithostatic_stress; % normal stress perpendicular to faults MPa
inputParameters.boundaryConditions.shear = shear_stress; % shear stress on faults MPa
% parameters tacked on to input file
fault_flaw = [];
fault_flaw.stiffS = shear_stiffness; % shear stiffness
fault_flaw.stiffN = normal_stiffness; % normal stiffness
fault_flaw.T = T; % Tensional strenght
fault_flaw.S_0 = S_0; % Internal shear strength
fault_flaw.cohes = sliding_cohesion; % sliding cohesion
fault_flaw.friction_s = static_friction; % static friction
fault_flaw.friction_d = dynamic_friction; % dynamic friction
fault_flaw.L = critical_slip_distance; % critical sliding distance
GROWInputParameters.fault_flaw = fault_flaw;
% software:
GROWInputParameters.GROW_perl_fileName = 'GROW_OCT_2017.pl';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ------------------- section functions ----------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check input
function [X,Y, plotOption] = checkinput(fileName, X,Y, input, numInput)
if isa(fileName, 'string')
error('input '' fileName must be of ype ''string''');
end
if numInput < 3
error('There must be at least three inputs: an X array, a Y array and a string specifying the file name')
end
if min(size(X)) ~=1
error('X array must be 1 by N or N by 1')
end
if min(size(Y)) ~=1
error('Y array must be 1 by N or N by 1')
end
if length(X) ~= length(Y)
error('X and Y arrays must be same length')
end
if length(X(1,:)) == length(X)
X = X';
end
if length(Y(1,:)) == length(Y)
Y = Y';
end
if ~isempty(input)
plotOption = input{1};
else
plotOption = 'long';
end
end
% check whether necessary files are in folder:
function check_directory(GROWInputParameters)
% check directory
checkFile(GROWInputParameters.GROW_perl_fileName);
checkFile('simple_shear_bc.pl');
checkFile('Wext.pl');
% compile fric2d if hasnt been done so already and move file into working
% directory
if (exist('fric2d', 'file') == 0)
% IMPORTANT: if fric
path2fric2dCompiler = 'GROW/fric2d_source_code';
disp('fric2d is not in the working directory, we will try to take care of it...')
currentFolder = pwd;
try %to compile fric 2d
cd(path2fric2dCompiler);
system('make');
cd(currentFolder);
system(['cp ', path2fric2dCompiler, '/fric2d fric2d']);
disp('got it')
catch
error(['the fric2d executable is not in working directory, a compiler was not found in the following directory: ', ...
path2fric2dCompiler]);
end
end
function checkFile(fileName)
if exist(fileName, 'file') ~= 2
errorMessage = [fileName, ' is not in directory, you need it to be'];
error(errorMessage)
end
end
end
%% 1) make the input file for the fric2D program
function [meanSegmentLength, X] = makeinputfile(fileName,X,Y, inputParameters)
% create profile:
profileName = [fileName,'_profile.txt']; % dont change this same name used in 6b
padding = 2;
boxAspectRatio = 0.3; % height to length ratio of the boundary box
[X, meanSegmentLength,boxSize] = makeprofile(profileName,X,Y,padding, boxAspectRatio);
% create boundary condition (creates a file named 'input')
% input ptspacing based on number of points
% simple_shear_bc.pl 0.06 0.15 0.0002
% (based on profile1)
input_section_fileName = 'input';
commandFormat = 'perl simple_shear_bc.pl %f %f %f %f %f %f %f %f';
command = sprintf(commandFormat , ...
boxSize(2) , ...
boxSize(1) , ...
meanSegmentLength , ...
inputParameters.material.E , ...
inputParameters.material.nu , ...
inputParameters.boundaryConditions.sigxx, ...
inputParameters.boundaryConditions.sigyy, ...
inputParameters.boundaryConditions.shear);
system(command);
% calculate number of endbit segments on either side of input fault
% geometry (magic number 4 is used to calulate half the distance from the
% tip of the fault geometry to the end side of the box)
numBufferSegments = floor(boxSize(1)/(padding*4*meanSegmentLength));
make_end_bits(X,Y, numBufferSegments, meanSegmentLength, boxSize, inputParameters);
% add tag line for the first fault
fault_header_fileName = 'fault_header.txt';
fid=fopen(fault_header_fileName,'w');
fprintf(fid, 'fault segment_0 yes no no\n');
fclose(fid);
% concatenate to the input file
command = sprintf('cat %s %s %s %s > %s', ...
input_section_fileName, ...
fault_header_fileName, ...
profileName, ...
'end_bits.txt', ...
[fileName, '.in'] );
system(command);
% more options...
% options that may need to be tweeked:
% observation points (change 'odata' and 'nolines')
% tolerance
% 1a) make the fault geometry based on series of x-y points, get the mean BEM
% size and determine the appropriate box size
function [newX, MEAN_SEGMENT_LENGTH,boxSize] = makeprofile(fileName,X,Y,padding, boxAspectRatio)
% prints profile defined by vertical vectors X and Y in line element form
% in the file specified by input fileName
% padding [dim1,dim2] defines the size of the boundary condition box as a function of
% the profile length and amplitude of the profile.
% *edits*:
% wed jul12 2017 changed the padding so that vertical box dimensions are
% instead a function of the profile length and determined by the
% 'boxAspectRatio'
RX = range(X);
RY = range(Y);
MEAN_SEGMENT_LENGTH = mean(sqrt((X(2:end)-X(1:(end-1))).^2 +(Y(2:end)-Y(1:(end-1))).^2));
MEAN_SEGMENT_LENGTH = round(MEAN_SEGMENT_LENGTH,2,'significant');
% boxSize = padding.*[RX,RY]; % not great because the box size changes as a
% function of topography, which is bad if set in a loop
% instead:
boxSize = zeros(1,2);
boxSize(1) = padding(1)*RX;
boxSize(2) = boxSize(1)*boxAspectRatio;
% round boxsize to nearest multiple of boundary element length
boxSize = ceil(boxSize/MEAN_SEGMENT_LENGTH)*MEAN_SEGMENT_LENGTH;
% errors/warnings related to profile in relation to box size
if max(abs(Y))>boxSize(2)*0.1
if max(abs(Y))>boxSize(2)/2
error('Profile topography is too large and intersect the boundary box, choose a larger aspect ratio for the box by changing the vatiable ''boxAspectRatio'' ')
end
warning('Profile topography is large, consider choosing a larger aspect ratio for the box by changing the vatiable ''boxAspectRatio'' ');
end
% round to increments of the length of elements
newX = X - min(X) + boxSize(1)/2 - RX/2; % shift X to be centered in the observation box
fid=fopen(fileName,'w');
for n = 1:(length(newX)-1)
fprintf(fid, ...
'1 %f %f %f %f 0 1.00E+10 0 0 0 0.6 0 0\n', ...
[newX(n), Y(n), newX(n+1), Y(n+1)]');
end
fprintf(fid, '\n');
fclose(fid);
end
% 1b) add end bits to input file according to numBuffer segments (figiting
% needed)
function [] = make_end_bits(X,Y,numBufferSegments, meanSegmentLength, BOXSIZE, inputParameters) % MAKE ThE FUNCTION CALL MATCh
intro = {'*FRACTURES'; ...
'*num xbeg ybeg xend yend kode bvs bvn bvs bvn'; ...
'*---- ---- ---- ---- ---- ---- ------ ---- ----- --------'};
tempFileName = 'end_bits.txt';
fileID = fopen(tempFileName,'w');
fprintf(fileID,'%s\r\n',intro{:});
fclose(fileID);
% find ends of the fault
X1 = X(1);
Y1 = Y(1);
Xend = X(end);
Yend = Y(end);
% find end of each end bit
start1 = [X1-numBufferSegments*meanSegmentLength,Y1];
end1 = [X1,Y1];
start2 = [Xend,Yend];
end2 = [Xend+numBufferSegments*meanSegmentLength,Yend];
modelHeight = BOXSIZE(2); % CHECK THIS
make_one_EB('left' ,tempFileName, start1,end1,numBufferSegments,meanSegmentLength, modelHeight, inputParameters);
make_one_EB('right',tempFileName, start2,end2,numBufferSegments,meanSegmentLength, modelHeight, inputParameters);
function [] = make_one_EB(EBName, fileName, EBStart, EBEnd,numBufferSegments, meanSegmentLength, MODELHEIGHT, inputParameters)
% function to make one end bit to avoid redundancy...
% make set of X coordinates for each end bit
% Left end bit (LEB):
EBX = linspace(EBStart(1),EBEnd(1),numBufferSegments+1);
% make start and end
EBXBeg = EBX(1:end-1);
EBXEnd = EBX(2:end);
EBY = linspace(EBStart(2),EBEnd(2),numBufferSegments+1);
% make start and end
EBYBeg = EBY(1:end-1);
EBYEnd = EBY(2:end);
% ... and midpoints
EBXMid = (EBXBeg+EBXEnd)/2;
EBYMid = (EBYBeg+EBYEnd)/2;
EBUs = calc_end_bit_displacement(EBXMid,EBYMid, ...
MODELHEIGHT, ...
inputParameters.material.E , ...
inputParameters.material.nu , ...
inputParameters.boundaryConditions.sigxx, ...
inputParameters.boundaryConditions.sigyy, ...
inputParameters.boundaryConditions.shear, ...
meanSegmentLength);
% header section
FID = fopen(fileName,'a');
fprintf(FID, 'fracture %s no no\n',EBName);
kode = 3;
bvs = 0;
bvn = 0;
% individual elements
for iElement = 1:numBufferSegments
ROW1 = [ 1, ...
EBXBeg(iElement), ...
EBYBeg(iElement), ...
EBXEnd(iElement), ...
EBYEnd(iElement), ...
kode, ...
bvs,...
bvn,...
EBUs(iElement),...
inputParameters.boundaryConditions.sigxx ];
fprintf(FID, '%f %f %f %f %f %f %f %f %f %f\n',ROW1);
end
fprintf(FID, '\n');
fclose(FID);
end
function [US] = calc_end_bit_displacement(x,y,h,E,nu,sigyy,sigxx,shear,elementSize)
% analytical solution for the displacement (u_s) on the end bits
% Input (1 by 1 or 1 by N)
% [x,y]: center coordinate of the endbit element
% E: Yougn's Mofulous (MPa) 1 by
% nu: Poisson's ratio
% [sigxx, sigyy]: Lithostatic stress (MPa)
% shear: Shear stress (MPa)
% elementSize: Size of boundary elements (2x slip) (m)
% analytical calculation
US = -(x./E) .* ((1-nu.^2).*-sigxx - nu.*(1+nu).*-sigyy) +...
( -shear .* (y+h/2) .*2*(1+nu)./E) - ...
elementSize/4;
end % function [US] = calc_end_bit_displacement(x,y,E,nu,sigyy,sigxx,shear,elementSize)
end % function [] = make_end_bits(X,Y,numBufferSegments, BVSTop, BVSBot)
end
%% 2) run fric2D
function runfric2D(fileName)
command = ['./fric2d -i ', ...
fileName,'.in -v -o ', fileName,'.out'];
system(command);
end
%% 3) extract output into a structure (uses getsection, assign2struct)
function [OUTPUT_STRUCT] = extractoutput(fileName)
% extracts numeric blocks of information out of output file and saves them
% into matlab structure with fields names corresponding (sometimes with
% small adjustments to header lines)
% input:
% fileName: name of the file in a string without the extension (here
% assumed to be .out). IMPORTANT: output file is assumed to have 9 data
% blocks. Namely 2 boudaryu line blocks, followed by 4 'end bit' blocks and
% 3 fault info blocks. The current state does not extract 'end bits'
% information.
% output:
% structure with with fields corresponding to blocks of data:
% 1) boundaryLines
% 2) boundaryConditions
% 3) faultinfo: information about fault geometry etc.
% 4) faultStress
% 5) faultStress2: additional fault stresses namely including the
% tangential stresses
% function calls:
% read_blocks_of_numeric_data
% assign2struct
% usage:
% extractoutput(sample_section_test)
% 2017-04-17: tested and functions properly using the file
% sample_section_test.out
outputFileName = [fileName,'.out'];
% load in block of data (TERRIBLE IDEA BUT FUCK IT):
expectedDataBlocks = 9; % this is here to make sure that I'm not fucking up the parsing of the output
dataBlocks = read_blocks_of_numerical_data(outputFileName, 10);
numBlocks = length(dataBlocks);
if numBlocks ~= expectedDataBlocks
error('the output is not of the expected size consider double-checking output parsing')
end
%extract the boundary lines:
boundaryLineBlock = dataBlocks(1);
boundaryLineHeaders = {'Elt' , ...
'xBeg' , ...
'yBeg' , ...
'xEnd' , ...
'yEnd' , ...
'length' , ...
'angle' , ...
'kode' , ...
'USstat' , ...
'UNstat' , ...
'USMonotonic' , ...
'UNMonotonic' };
boundaryLines = assign2struct(boundaryLineBlock{:},boundaryLineHeaders);
% extract boundary conditions:
boundaryConditionBlock = dataBlocks(2);
boundaryConditionHeaders = { 'Elt' , ...
'DS' , ...
'USminus' , ...
'USplus' , ...
'DN' , ...
'UNminus' , ...
'UNplus' , ...
'SigmaS' , ...
'SigmaN' };
boundaryConditions = assign2struct(boundaryConditionBlock{:},boundaryConditionHeaders);
% exctract fault information
faultInfoBlock = dataBlocks(7);
faultInfoHeaders = boundaryLineHeaders;
faultInfo = assign2struct(faultInfoBlock{:}, faultInfoHeaders);
% exctract fault information: stress conditions
faultStressBlock = dataBlocks(8);
faultStressHeaders = boundaryConditionHeaders;
faultStress = assign2struct(faultStressBlock{:}, faultStressHeaders);
% exctract fault information: stress conditions
faultStress2Block = dataBlocks(9);
faultStress2Headers = { 'BE' , ...
'Tangential' , ...
'maxPrinctop' , ...
'angleTop' , ...
'tangential' , ...
'maxPrincBot' , ...
'angleBot' };
faultStress2 = assign2struct(faultStress2Block{:}, faultStress2Headers);
% extract endbit geometry:
leftEndBitBlock = dataBlocks(3);
rightEndBitBlock = dataBlocks(5);
endBitHeaders = boundaryLineHeaders;
leftEndBit = assign2struct(leftEndBitBlock{:}, endBitHeaders);
rightEndBit = assign2struct(rightEndBitBlock{:}, endBitHeaders);
OUTPUT_STRUCT = [];
OUTPUT_STRUCT.boundaryLines = boundaryLines;
OUTPUT_STRUCT.boundaryConditions= boundaryConditions;
OUTPUT_STRUCT.faultInfo = faultInfo;
OUTPUT_STRUCT.faultStress = faultStress;
OUTPUT_STRUCT.faultStress2 = faultStress2;
OUTPUT_STRUCT.leftEndBit = leftEndBit;
OUTPUT_STRUCT.rightEndBit = rightEndBit;
end
%% 4) assess failure based on model output
function [shearFailureCoord,tensileFailureCoord, ...
coulombCriterion, tensileCriterion] = assessfailure(X, Y, OUTPUT_DATA_STRUCT, ...
inputParameters)
% assessed the failure of BE's along the fault profile. This is done by 1)
% using the structure output produced by the extractouput function to
% determine the stress tensor, 2) finding the prinple stesses (the
% eigenvalues of the stress tensor) and 3) assessing both shear
% (mohr-coulomb) and tensile failure criterions. 4) finding the
% corresponding locations on the profile
% INPUT:
% OUTPUT_DATA_STRUCT: output data extracted using extractoutput
% OUTPUT:
% shearFailureCoord: 2 by N array with x-y coordinates of elements that failed in shear
% tensileFailureCoord: 2 by N array with x-y coordinates of elements that failed in tnesion
% coulombCriterion: 1 by N array coulomb stress
% (tau_m-[\sigma_m*sin(\phi) + c*cos(\phi)])
% where
% \sigma_m = (\sigma_1+\sigma_3)/2
% \tau_m = (\sigma_1-\sigma_3)/2
% tensileCriterion: 1 by N array tensile criterion
% (-\sigma_3-c)
% if the latter 2 are positive, the fault elements are expected to fail
% function calls:
% gettensor
% usage:
% assessfailure(outputDataStruct)
% 2017-04-17: tested and functions properly using the file
% sample_section_test.out processed with function extractoutput (did not
% test the physical viability of results)
c = inputParameters.material.c; % cohesion
phi = inputParameters.material.phi; % angle of internal friction
[sigma11, sigma22, sigma12] = gettensor(OUTPUT_DATA_STRUCT,'top');
[shearFailureCoordTop, tensileFailureCoordTop] = find_failure_points;
[sigma11, sigma22, sigma12] = gettensor(OUTPUT_DATA_STRUCT,'bottom');
[shearFailureCoordBottom, tensileFailureCoordBottom] = find_failure_points;
shearFailureCoord = [shearFailureCoordTop; shearFailureCoordBottom];
tensileFailureCoord = [tensileFailureCoordTop; tensileFailureCoordBottom];
function [shearFailureCoord, tensileFailureCoord] = find_failure_points()
% returns sorted (as a function of magnitude) coordinates for failure
% both in shear and in tensile stress)
% uses variables: sigma11, sigma22, sigma12, OUTPUT_DATA_STRUCT
numBE = length(sigma11);
coulombCriterion= zeros(numBE,1);
tensileCriterion= zeros(numBE,1);
for iBE = 1:numBE
tensor = [sigma11(iBE), sigma12(iBE); ...
sigma12(iBE), sigma22(iBE)];
eigenValues = eig(tensor);
% DOING GEOLOGY FOR A MOMENT HERE
sigma1 = -eigenValues(2);
sigma3 = -eigenValues(1);
tau_m = (sigma1-sigma3)/2;
sigma_m = (sigma1+sigma3)/2;
coulombCriterion(iBE) = -tau_m + (sigma_m*sind(phi)+c*cosd(phi));