data/tools/support/maximum.py~

+import numpy as np
+
+def maximum(x) :
+    size = x.shape 
+    res = np.zeros(size[0])
+
+    for time_step in range(size[0]) :
+        for vertex in range(size[1]) :
+            res[time_step] = max(res[time_step], max(x[time_step, vertex]))
+
+    return res
\ No newline at end of file

data/tools/support/norm.py

+import numpy as np
+
+
+def norm(x):
+    size = x.shape
+    res = np.zeros((size[0], size[1]))
+
+    for time_step in range(size[0]):
+        for vertex in range(size[1]):
+            res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])
+
+    return res

data/tools/support/norm.py~

+import numpy as np
+
+def norm(x) :
+    size = x.shape 
+    res = np.zeros(size[0], size[1])
+
+    for time_step in range(size[0]) :
+        for vertex in range(size[1]) :
+            res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])
+
+    return res 

data/tools/support/slip_beginnings.py

+import array as ar
+
+def slip_beginnings(x):
+    # returns indicies i for which x[i-1]=0 and x[i]>0
+    starts = ar.array('i', (0 for i in range(len(x))))
+    prev = False
+
+    length = 0
+    for i in range(len(x)):
+        if (not prev and x[i]):
+            starts[length] = i
+            length += 1
+        prev = x[i]
+    return starts[:length]

data/tools/support/slip_beginnings.py~

+import numpy as np
+
+def slip_beginnings(x) :
+    # returns indicies i for which x[i-1]=0 and x[i]>0  
+    starts = []
+    prev = False
+    for i in range(len(x)) :
+        if (not prev and x[i]) :
+            starts.append(i)
+            prev = x[i]
+    return starts

data/tools/support/slip_endings.py

+import array as ar
+
+def slip_endings(x):
+    # returns indicies i for which x[i-1]>0 and x[i]=0
+    ends = ar.array('i', (0 for i in range(len(x))))
+    prev = False
+
+    length = 0
+    for i in range(len(x)):
+        if (prev and not x[i]):
+            ends[length] = i
+            length += 1
+        prev = x[i]
+    return ends[:length]
\ No newline at end of file

data/tools/support/slip_endings.py~

+import numpy as np
+
+def slip_endings(x) :
+    # returns indicies i for which x[i-1]>0 and x[i]=0
+    starts = np.array([])
+    prev = False
+    for i in range(len(x)) :
+        if (prev and not x[i]) :
+            starts.append(i)
+            prev = x[i]
+    return starts
\ No newline at end of file

data/tools/support/slip_rates.py

+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+
+SORT_COORD = 0
+
+def local_to_global(coords, SORT_COORD=0):
+    return coords[:, SORT_COORD].argsort()
+
+def slip_rates(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIAL_COORDS=1):
+    body = 'body' + str(body_ID)   # 'frictionalBoundary'     'body' + str(body_ID)
+
+    coords = np.array(h5file[body + '/coordinates'])
+    if len(patch) == 0:
+        patch = np.linspace(0, len(coords)-1, len(coords), endpoint=True, dtype='int8')
+
+    # local to global
+    local_global = local_to_global(coords, SORT_COORD) 
+    local_global = local_global[patch]
+
+    patch_size = coords[local_global[-1], SORT_COORD] - coords[local_global[0], SORT_COORD]
+
+    # read time
+    time = np.array(h5file['relativeTime']) * FINAL_TIME
+    time = np.delete(time, 0)
+    if len(interval) == 0:
+        interval = [0, FINAL_TIME]
+    t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
+    time = time[t]
+
+    fig = plt.figure()
+
+    # velocity data
+    v = abs(np.array(h5file[body + '/velocity']))
+    v_t = v[t,:,TANGENTIAL_COORDS]
+    v_tx = v_t[:,local_global]
+
+    # plot
+    ax_slip = fig.add_subplot(1, 1, 1)
+    pos = ax_slip.imshow(v_tx, interpolation='nearest', cmap=cm.Greys_r, aspect='auto', origin='lower', extent=[0, patch_size, time[0], time[-1]])
+    #ax_slip.set_ylim([time[-1], time[0]])
+    ax_slip.set_ylabel('time [s]')
+    ax_slip.set_xlabel('horizontal coordinate [m]')
+    fig.colorbar(pos, ax=ax_slip)
+
+    #ax_slip.set_yscale('log')
+   # ax_slip.set_ylim([1e-6,1e-2])
+    #-------------------------
+
+    fig.canvas.draw()
\ No newline at end of file

data/tools/support/trapezoidal.py

+# x: vector of domain values
+# y: vector of function values f(x)
+def trapezoidal(x, y) :
+    # returns integral of (x, y = f(x)) calculated with trapezoidal rule 
+	res = 0.0
+	for i in range(1, len(x)) :
+		res += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2.0
+	return res
\ No newline at end of file

data/tools/support/trapezoidal.py~

+# x: vector of domain values
+# y: vector of function values f(x)
+def trapezoidal(x, y) :
+	# returns integral of (x, y = f(x)) calculated with trapezoidal rule 
+	res = 0.0
+	for i in range(1, len(x)) :
+		res += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2.0
+	return res
\ No newline at end of file

data/tools/support/write_contours.py

+writeContours <- function (contours, level, file) {
+  file.create(file)
+  for (cl in contours[sapply(contours, function(x) x$level) == level]) {
+    conn <- file(file, 'a')
+    write.table(cbind(cl$x, cl$y, level),
+                file = file, row.names = FALSE, col.names = FALSE,
+                append = TRUE)
+    writeLines("\n", conn)
+    close(conn)
+  }
+}

data/tools/truncation.py

+import configparser as cp
+import os
+import numpy as np
+import csv
+import h5py
+import matplotlib.pyplot as plt
+
+from debug.outliers import outliers
+from debug.friction import truncated_friction
+from debug.state import aging_law
+from debug.diffplot import diffplot
+
+from support.maximum import maximum
+from support.norm import norm
+from support.find_quakes import find_quakes
+from support.slip_beginnings import slip_beginnings
+from support.slip_endings import slip_endings
+from support.max_distance import max_distance
+
+from support.io import read_h5file
+from support.io import read_params
+
+from support.iterations import iterations
+from support.friction_stats import friction_stats
+
+def build_time(h5file, final_time, interval = []):
+    # read time
+    time = np.array(h5file['relativeTime']) * FINAL_TIME
+    time = np.delete(time, 0)
+    
+    if len(interval) == 0:
+        interval = [0, FINAL_TIME]
+    
+    t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
+    return t, time[t]
+
+def build_patch(coords, percentage = 1.0):
+    x_coords = coords[:, 0]
+    xmin = np.min(x_coords)
+    xmax = np.max(x_coords)
+    delta_x = (1 - percentage)*(xmax - xmin)/2
+
+    xmin = xmin + delta_x
+    xmax = xmax - delta_x
+
+    return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]
+
+# read h5 output file
+h5file = read_h5file()
+print(list(h5file.keys()))
+
+# read problem parameters
+params = read_params('foam.cfg')
+print(params)
+
+TANGENTIAL_COORDS = 1
+FINAL_TIME = params['finalTime']
+NBODIES = params['bodyCount']
+
+t, time = build_time(h5file, FINAL_TIME, [0, FINAL_TIME])
+
+for body_ID in range(NBODIES):
+    body = 'body' + str(body_ID)
+
+    if body not in h5file:
+        continue
+
+    coords = np.array(h5file[body + '/coordinates'])
+    patch = build_patch(coords, 1.0)
+
+    fig = plt.figure()
+
+    # velocity and state data
+    v = abs(np.array(h5file[body + '/velocity']))[:,:,TANGENTIAL_COORDS]
+    alpha = np.array(h5file[body + '/state'])
+
+    vmin = 0 #params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )
+    
+    v_truncation = (v<=vmin)[t,:]
+    v_truncated = np.sum(v_truncation, axis=1)
+
+    v_comp = (np.abs(v-vmin)<1e-14)[t,:]
+    v_small = np.sum(v_comp, axis=1)
+
+    weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
+    second_deriv_large = (v<0)[t,:] #((- weighted_normal_stress[None, ...] * params['a'] / v)>1e10)[0,t,:]
+    second_deriv = np.sum(second_deriv_large, axis=1) 
+
+    total_truncated = np.sum(v_truncation | v_comp | second_deriv_large, axis=1)
+
+    # plot v_small
+    ax_v_small = fig.add_subplot(4, 1, 1)
+    ax_v_small.plot(time, v_small, color='black', linestyle='-')
+    ax_v_small.set_ylabel('# v-vmin < 1e-14')
+    #-------------------------
+
+    # plot v_truncated
+    ax_v_truncated = fig.add_subplot(4, 1, 2  )
+    ax_v_truncated.plot(time, v_truncated, color='black', linestyle='-')
+    ax_v_truncated.set_ylabel('# v <= vmin')
+    #-------------------------
+
+    # plot second_deriv_large
+    ax_second_deriv = fig.add_subplot(4, 1, 3)
+    ax_second_deriv.plot(time, second_deriv, color='black', linestyle='-')
+    ax_second_deriv.set_ylabel('# second_deriv large')
+    #-------------------------
+
+    # plot total_truncated
+    ax_truncated = fig.add_subplot(4, 1, 4)
+    ax_truncated.plot(time, total_truncated, color='black', linestyle='-')
+    ax_truncated.set_ylabel('# total truncated')
+    #-------------------------
+plt.show()
+
+h5file.close()
\ No newline at end of file

debugging.m

+% householder reflection
+
+v0 = [1; 0];
+v0 = v0./norm(v0);
+
+v1 = [1; -1];
+v1 = v1./norm(v1);
+
+v = 1/2.*(v1-v0);
+vnorm = norm(v);
+
+H = eye(2) - 2/(vnorm^2).*(v*v')
+
+A = [...
+];
+
+f = [...
+0;
+-4.33681e-19;
+0;
+1.0842e-19;
+-2.66671e-06;
+-4.94914e-07;
+0;
+-2.71051e-20;
+-4.19923e-06;
+-4.21374e-07;
+-3.92235e-06;
+-8.64817e-07;
+0;
+-5.42101e-20;
+-4.96301e-06;
+-4.13502e-06;
+0;
+0;
+-7.03515e-06;
+-5.84881e-06;
+-4.93309e-06;
+-2.37804e-06;
+0;
+0;
+-2.90755e-06;
+-1.23334e-06;
+0;
+0;
+-2.2449e-06;
+-8.21563e-07;
+0;
+0;
+0;
+0;
+0;
+0;
+-7.24652e-06;
+-8.25293e-07;
+0;
+2.4104e-06;
+-8.19294e-06;
+-1.18617e-06;
+-5.58223e-06;
+-1.69348e-06;
+0;
+1.88001e-06;
+-1.03387e-05;
+-8.22886e-06;
+0;
+0;
+-1.07258e-05;
+-8.62344e-06;
+-9.76077e-06;
+-4.93954e-06;
+-9.10261e-06;
+-7.35307e-06;
+-4.01973e-06;
+-2.38351e-06;
+-5.92659e-06;
+-4.23858e-07;
+-7.9911e-06;
+-4.04139e-06;
+-3.80594e-06;
+-3.06903e-06;
+-3.92924e-06;
+-2.71681e-06;
+-5.12975e-06;
+-1.37642e-06;
+-4.79892e-06;
+-1.14935e-06;
+-7.19726e-06;
+-3.54967e-06;
+-6.62755e-06;
+-6.1499e-07;
+0;
+1.28701e-05;
+0;
+0;
+-1.09944e-05;
+-1.92221e-06;
+0;
+1.17853e-05;
+-1.21456e-05;
+-2.20784e-06;
+-8.3682e-06;
+-2.01449e-06;
+0;
+0;
+-1.43122e-05;
+-1.00938e-05;
+0;
+0;
+-1.48374e-05;
+-1.03947e-05;
+-1.3684e-05;
+-6.31825e-06;
+-1.35318e-05;
+-9.3703e-06;
+-5.76679e-06;
+-1.88878e-06;
+-9.18524e-06;
+-1.51544e-06;
+-1.16296e-05;
+-5.64398e-06;
+-5.20301e-06;
+-2.13559e-06;
+-5.55432e-06;
+-2.017e-06;
+-7.96841e-06;
+-1.7895e-06;
+-7.95375e-06;
+-1.76175e-06;
+-1.0596e-05;
+-5.24553e-06;
+-1.02345e-05;
+-1.74422e-06;
+-2.16064e-05;
+-1.85059e-06;
+-2.24752e-05;
+-3.34957e-07;
+-1.34774e-05;
+-2.51353e-06;
+-2.18312e-05;
+-1.04202e-06;
+-1.99898e-05;
+-1.463e-06;
+-1.94256e-05;
+-2.11821e-06;
+-5e-05;
+5.03269e-07;
+-2.00631e-05;
+-2.41829e-06;
+-2.11637e-05;
+-1.12298e-06;
+-1.61508e-05;
+-1.11752e-05;
+-1.7267e-05;
+-6.9337e-06;
+-1.43931e-05;
+-6.46903e-06;
+-1.70354e-05;
+-2.45177e-06;
+-5e-05;
+-2.50465e-06;
+-5e-05;
+-1.21588e-05;
+-2.01431e-05;
+-1.19797e-05;
+];
+
+d = 0.00075745;
+A2 = [...
+];
+
+f2 = [...
+0;
+0;
+0;
+0;
+-3.01994e-06;
+-1.76298e-07;
+0;
+0;
+-4.78204e-06;
+-2.0995e-07;
+-5.01951e-06;
+1.22653e-06;
+0;
+0;
+-4.98476e-06;
+-3.57931e-06;
+0;
+0;
+-7.06572e-06;
+-4.97786e-06;
+-5.40025e-06;
+-1.91591e-06;
+0;
+0;
+-3.26344e-06;
+-9.4196e-07;
+0;
+0;
+-3.47107e-06;
+1.69693e-06;
+0;
+0;
+0;
+0;
+0;
+0;
+-8.36624e-06;
+-3.21786e-07;
+0;
+0;
+-9.19771e-06;
+-4.8646e-07;
+-9.46319e-06;
+2.13218e-06;
+0;
+0;
+-1.03261e-05;
+-6.78209e-06;
+0;
+0;
+-1.09634e-05;
+-7.09926e-06;
+-1.01615e-05;
+-3.77678e-06;
+-9.09151e-06;
+-6.13484e-06;
+-9.06507e-06;
+4.28268e-06;
+-6.83374e-06;
+-2.4706e-07;
+-8.51881e-06;
+-3.22169e-06;
+-7.83369e-06;
+3.49294e-06;
+-8.45665e-06;
+3.90404e-06;
+-7.913e-06;
+1.91502e-06;
+-7.13457e-06;
+1.77604e-06;
+-7.71821e-06;
+-2.87949e-06;
+-7.58455e-06;
+-3.19117e-07;
+0;
+0;
+0;
+0;
+-1.17905e-05;
+-6.57623e-07;
+0;
+0;
+-1.26738e-05;
+-7.75536e-07;
+-1.28663e-05;
+2.59111e-06;
+0;
+0;
+-1.33726e-05;
+-8.23012e-06;
+0;
+0;
+-1.40382e-05;
+-8.4827e-06;
+-1.3445e-05;
+-4.60248e-06;
+-1.21277e-05;
+-7.7075e-06;
+-1.20933e-05;
+5.59837e-06;
+-1.00511e-05;
+-5.83173e-07;
+-1.16989e-05;
+-4.19969e-06;
+-1.07938e-05;
+5.03904e-06;
+-1.14551e-05;
+5.33497e-06;
+-1.1057e-05;
+2.39173e-06;
+-1.01977e-05;
+2.22842e-06;
+-1.08385e-05;
+-3.93909e-06;
+-1.09119e-05;
+-6.59537e-07;
+-1.51572e-05;
+6.80143e-06;
+-1.97651e-05;
+7.92292e-06;
+-1.3582e-05;
+-8.81729e-07;
+-1.73372e-05;
+7.52177e-06;
+-1.63794e-05;
+2.79237e-06;
+-1.33659e-05;
+2.38778e-06;
+-5e-05;
+7.93536e-06;
+-1.92072e-05;
+-9.96979e-07;
+-1.93357e-05;
+3.15138e-06;
+-1.54595e-05;
+-9.11795e-06;
+-1.68142e-05;
+-5.1397e-06;
+-1.4113e-05;
+-4.66534e-06;
+-1.63265e-05;
+-9.49749e-07;
+-5e-05;
+-1.0631e-06;
+-5e-05;
+-1.02021e-05;
+-1.98778e-05;
+-9.99727e-06;
+];
+
+d2 = 0;
+
+sum(sum(abs(f-f2)>1e-14))
+
+
+sum(sum(abs(A-A2)>1e-14))
+
+
+alpha0 = 100;
+V = 1e-5; %1e-20;
+V0 = 5e-5;
+L = 2.25e-5;
+
+alpha = @(t) log(exp(alpha0 - V.*t./L ) - V0 * (exp(-t./L.*V)-1)./V);
+
+x = [0:0.001:0.1];
+alpha(x)
+plot(x, alpha(x));
+
+%abs(d-d2)
+
+% x = [...
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% -1.2337e-10;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% 0;
+% -1.2337e-10;
+% 0;
+% -1.2337e-10;
+% 0;
+% 0;
+% 0;
+% ];
+% 
+% ignore = [...
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0; 
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 1;
+% 1; 
+% 0;
+% 0; 
+% 1; 
+% 1; 
+% 0;
+% 0; 
+% 1;
+% 1; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 1;
+% 0; 
+% 0;
+% 0; 
+% 0;
+% 0; 
+% 0; 
+% 0; 
+% 0;
+% 0; 
+% 0; 
+% 0; 
+% 0; 
+% 0; 
+% 1;
+% 0; 
+% 1;
+% 0; 
+% 0;
+% 0; 
+% ];
+% 
+% newA = A;
+% newf = f;
+% for i=1:length(ignore)
+%     if (ignore(i)==1)
+%         for j=1:length(ignore)
+%             newA(i,j) = (i==j);
+%         end
+%         newf(i) = x(i);
+%     end
+% end
+% newA
+% newf
+
+% ignore = zeros(1, length(f));
+% for i=1:length(ignore)
+%     if (A(i,i)==0)
+%         ignore(1, i) = 1;
+%     end
+% end
+% 
+% indices = [];
+% for i=1:length(ignore)
+%     if (ignore(i)==0)
+%         indices = [indices i];
+%     end
+% end
+% 
+% len = length(indices);
+% newA = zeros(len, len);
+% newf = zeros(1, len);
+% for i=1:len
+%     newf(i) = f(indices(i));
+% 
+%     for j=1:len
+%         newA(i,j) = A(indices(i), indices(j));
+%     end
+% end
+% newf;
+% newA
+
+%inv(A);
+
+% sol = newA\newf
+% 
+% tnnmgSol = [...
+%     -1.03349;
+% -0.0502635;
+% 0.0161336;
+% -0.00852914;
+% -0.0130659;
+% 0.00660115;
+% -0.434915;
+% -0.0128241;
+% -0.00893317;
+% 0.00884256;
+% -0.0110232;
+% -0.10936;
+% 1.14978;
+% 0.0581132;
+% 0.019598;
+% 0.190758;
+% 0.478393;
+% -0.00433049;
+% 0.0513446;
+% 0.353218;
+% -0.00482644;
+% 0.11758;
+% 0;
+% 0;
+% -0.013677;
+% 0.0512563;
+% 0;
+% 0;
+% -0.00514197;
+% -0.117664;
+% 0;
+% 0;
+% 0.130034;
+% 0.574041;
+% 0.0136965;
+% 0.214953;
+% 0.0335986;
+% -0.00703477;
+% 0.0438202;
+% 0.346506;
+% -0.000896805;
+% 0.127837;
+% 0.041946;
+% 0.217409;
+% -1.2337e-10;
+% 0.112528;
+% -0.00828618;
+% 0.00453345;
+% -0.00898709;
+% 0.0781683;
+% 0.092355;
+% -0.556143;
+% -0.00443686;
+% -0.140953;
+% 0.0424371;
+% -0.250082;
+% -0.00336878;
+% 0.00105471;
+% -1.2337e-10;
+% 0.0188452;
+% -1.2337e-10;
+% -0.14111;
+% 0.00298499;
+% -0.190904;
+% ];
+% 
+% sol-tnnmgSol
\ No newline at end of file

dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -5,4 +5,4 @@
 Module: dune-tectonic
 Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Depends: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/CMakeLists.txt

+add_subdirectory("data-structures")
+add_subdirectory("factories")
+add_subdirectory("io")
+add_subdirectory("multi-threading")
+add_subdirectory("problem-data")
+add_subdirectory("spatial-solving")
+add_subdirectory("tests")
+add_subdirectory("time-stepping")
+add_subdirectory("utils")
+
+add_custom_target(tectonic_dune SOURCES
+  assemblers.hh
+  assemblers.cc 
+  explicitgrid.hh
+  explicitvectors.hh
+  gridselector.hh
+)
+
+#install headers
 install(FILES
-  body.hh
-  frictiondata.hh
-  frictionpotential.hh
-  globalfrictiondata.hh
-  globalfriction.hh
-  globalratestatefriction.hh
-  gravity.hh
-  localfriction.hh
-  minimisation.hh
-  myblockproblem.hh
-  mydirectionalconvexfunction.hh
-  quadraticenergy.hh
-  tectonic.hh
-  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
+  assemblers.hh
+  explicitgrid.hh
+  explicitvectors.hh
+  gridselector.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)