data/tools/2d-dip-contours.R

+source('tools/support/findQuakes.R')
+source('tools/support/writeContours.R')
+Rcpp::sourceCpp('tools/support/trapezoidal.cpp')
+
+finalTime              <- 1000 # s
+specialTrenchDistances <- c(0.15,0.30,0.45) # m
+convergenceVelocity    <- 5e-5 # m/s
+
+last      <- function(x) x[[length(x)]]
+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (basedir in c("rfpitol=100e-7")) {
+  dir          <- file.path(directories[['simulation']],
+                            '2d-lab-fpi-tolerance', basedir)
+  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+  relativeTime <- h5file['relativeTime'][]
+  realTime     <- finalTime * relativeTime
+
+  calcMask     <- relativeTime >= 0.7
+  calcTime     <- realTime[calcMask]
+  calcDuration <- last(calcTime) - calcTime[[1]]
+  calcRange    <- range(which(calcMask))
+  calcLength   <- sum(calcMask)
+
+  ## We are interested in an enlarged time range around actual events here,
+  ## (and no other quantities!) hence we pass a very low velocity here.
+  quakes           <- findQuakes(1e-6 + convergenceVelocity,
+                                 h5file['/frictionalBoundary/velocity'][,,],
+                                 indices = calcRange[1]:calcRange[2])
+  quakes$beginning <- realTime[quakes$beginningIndex]
+  quakes$ending    <- realTime[quakes$endingIndex]
+  quakes$duration  <- quakes$ending - quakes$beginning
+  numQuakes        <- nrow(quakes)
+
+  relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
+                               quakes[[numQuakes,  'ending']]), f=0.02)
+  threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+
+  plotMask   <- threeQuakeTimeMask
+  plotTime   <- realTime[plotMask]
+  plotRange  <- range(which(plotMask))
+  plotLength <- sum(plotMask)
+
+  write(relaxedTime[[1]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'min', 'threequakes', basedir),
+                         'tex')))
+  write(relaxedTime[[2]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'max', 'threequakes', basedir),
+                         'tex')))
+
+  surfaceCoordinates    <- h5file['/surface/coordinates'][]
+  surfaceTrenchDistance <- surfaceCoordinates[,1] / cos(atan(.27))
+  perm                  <- order(surfaceTrenchDistance)
+
+  displacement <- h5file['/surface/displacement']
+  # This is the displacement imposed through the Dirichlet condition at the last node
+  displacementOffset <- displacement[calcRange[1]:calcRange[2],last(perm),1]
+  balancedSurfaceDisplacement <- matrix(nrow = plotLength, ncol = displacement@dim[2])
+  for (k in 1:displacement@dim[2]) {
+    d <- displacement[calcRange[1]:calcRange[2],k,1:2]
+    d[,,1] <- d[,,1] - displacementOffset
+    # We're in a tilted coordinate system
+    dv <- -sin(atan(.27))*d[,,1] + cos(atan(.27))*d[,,2]
+    meanVertSurfaceDisplacement  <- trapezoidal(calcTime, dv) / calcDuration
+    balancedSurfaceDisplacement[,k] <- (dv - meanVertSurfaceDisplacement)[
+      (plotRange[1] - calcRange[1] + 1):(plotRange[2] - calcRange[1] + 1)
+    ]
+  }
+
+  ## Interpolate velocity to special points from surrounding velocities
+  {
+    data <- matrix(nrow=plotLength, ncol=1+length(specialTrenchDistances))
+    data[,1] <- plotTime
+    for (k in seq(plotLength)) {
+      interpolant <- approxfun(surfaceTrenchDistance[perm],
+                               balancedSurfaceDisplacement[k,perm])
+      for (i in seq(specialTrenchDistances)) {
+        specialTrenchDistance <- specialTrenchDistances[[i]]
+        data[k,i+1] <- interpolant(specialTrenchDistance)
+      }
+    }
+    singleDataName <- file.path(directories[['output']],
+                                paste.(pasteColon('dip-single-points', basedir),
+                                       'csv'))
+    write.csv(data, singleDataName, row.names = FALSE)
+  }
+  h5::h5close(h5file)
+
+  printlevels <- c('-40','-20','-10','-5','-2.5','0','2.5','5','10','20','40')
+  levels      <- as.numeric(printlevels) * 1e-6
+
+  ret <- contourLines(plotTime,
+                      surfaceTrenchDistance[perm],
+                      balancedSurfaceDisplacement[,perm],
+                      levels=levels)
+
+  for (i in seq(printlevels))
+    writeContours(ret, levels[[i]],
+                  file.path(directories[['output']],
+                            paste.(pasteColon('2d-dip-contours', basedir,
+                                              'level', printlevels[[i]]),
+                                   'tex')))
+}

data/tools/2d-event-writer.R

+source('tools/support/findQuakes.R')
+
+finalTime              <- 1000 # s
+convergenceVelocity    <- 5e-5 # m/s
+criticalVelocity       <- 1000e-6 + convergenceVelocity
+
+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (basedir in c("rfpitol=100e-7")) {
+  dir          <- file.path(directories[['simulation']],
+                            '2d-lab-fpi-tolerance', basedir)
+  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+  relativeTime <- h5file['relativeTime'][]
+  realTime     <- finalTime * relativeTime
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  quakes <- findQuakes(criticalVelocity, velocityProxy,
+                       indices = 1:dim(velocityProxy)[1], 1)
+
+  basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]
+
+  maximumVelocities<- maxVelocity(velocityProxy[,,1])
+  displacement     <- h5file['/frictionalBoundary/displacement']
+  numVertices      <- dim(displacement)[[2]]
+  for (quakeID in seq(nrow(quakes))) {
+    qb <- quakes[[quakeID,'beginningIndex']]
+    qe <- quakes[[quakeID,'endingIndex']]
+    localSlippingTimes <- abs(velocityProxy[qb:qe,,1][,,1]) > criticalVelocity
+
+    qd   <- displacement[qb:qe,,1]
+    slip <- vector(mode='numeric', length=numVertices)
+    for (i in seq(numVertices)) {
+      if (any(localSlippingTimes[,i])) {
+        begs <- positiveStarts(localSlippingTimes[,i])
+        ends <- negativeStarts(localSlippingTimes[,i])
+        slip[i]<- sum(qd[ends,i,] - qd[begs,i,])
+      }
+    }
+    quakes[quakeID,'peakSlip']     <- max(abs(slip))
+    quakes[quakeID,'peakSlipRate'] <- max(maximumVelocities[qb:qe])
+
+    maxRuptureWidth <- 0
+    for (ts in seq(dim(localSlippingTimes)[[1]])) {
+      st <- localSlippingTimes[ts,]
+      if (!any(st))
+        next;
+
+      slippingCoordinates <- basalCoordinates[st,1] # x-coordinate only
+      maxRuptureWidth <- max(maxRuptureWidth, diff(range(slippingCoordinates)))
+    }
+    quakes[quakeID,'ruptureWidth'] <- maxRuptureWidth
+  }
+  h5::h5close(h5file)
+
+  quakes$beginning <- realTime[quakes$beginningIndex]
+  quakes$ending    <- realTime[quakes$endingIndex]
+
+  write.csv(quakes[c('beginning','ending',
+                     'peakSlip','peakSlipRate',
+                     'ruptureWidth')],
+            row.names = FALSE, quote = FALSE,
+            file = file.path(directories[['output']],
+                             paste.(pasteColon('events', basedir), 'csv')))
+}

data/tools/2d-fpi-tolerance.R

+rfpitols <- c('1e-7', '2e-7', '3e-7', '5e-7',
+              '10e-7', '20e-7', '30e-7', '50e-7',
+              '100e-7', '200e-7', '300e-7', '500e-7',
+              '1000e-7', '2000e-7', '3000e-7', '5000e-7',
+              '10000e-7', '20000e-7', '30000e-7', '50000e-7',
+              '100000e-7')
+numEntries <- length(rfpitols)
+
+data <- data.frame(row.names = rfpitols,
+                   tol = rep(NA, numEntries),
+                   time = rep(NA, numEntries),
+                   fpi = rep(NA, numEntries),
+                   mg = rep(NA, numEntries))
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (rfpitol in rfpitols) {
+  basedir <- paste('rfpitol', rfpitol, sep='=')
+  dir     <- file.path(directories[['simulation']],
+                       '2d-lab-fpi-tolerance', basedir)
+  h5file  <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+
+  data[rfpitol,'tol'] <- as.numeric(rfpitol)
+
+  relativeTimeProxy <- h5file['/relativeTime']
+  relativeTimeLen <- relativeTimeProxy@dim
+  data[rfpitol,'time'] <- relativeTimeProxy[relativeTimeProxy@dim]
+
+  ## FIXME: why do we drop the first entry?
+  fixedPointIterationsProxy <- h5file["/iterations/fixedPoint/total"]
+  fixedPointIterationsLen <- fixedPointIterationsProxy@dim
+  data[rfpitol,'fpi'] <- sum(fixedPointIterationsProxy[2:fixedPointIterationsLen])
+
+  multiGridIterationsProxy <- h5file["/iterations/multiGrid/total"]
+  multiGridIterationsLen <- multiGridIterationsProxy@dim
+  data[rfpitol,'mg'] <- sum(multiGridIterationsProxy[2:multiGridIterationsLen])
+  h5::h5close(h5file)
+}
+
+write.csv(data, file.path(directories[['output']], 'fpi-data.csv'),
+          row.names = FALSE, quote = FALSE)

data/tools/2d-performance.R

+source('tools/support/findQuakes.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+
+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (basedir in c("rfpitol=100e-7")) {
+  dir          <- file.path(directories[['simulation']],
+                            '2d-lab-fpi-tolerance', basedir)
+  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+  relativeTime <- h5file['relativeTime'][]
+  realTime     <- finalTime * relativeTime
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  ## We are interested in an enlarged time range around actual events here,
+  ## (and no other quantities!) hence we pass a very low velocity here.
+  quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
+                       indices = 1:dim(velocityProxy)[1], 1)
+  quakes$beginning <- realTime[quakes$beginningIndex]
+  quakes$ending    <- realTime[quakes$endingIndex]
+  quakes$duration  <- quakes$ending - quakes$beginning
+  numQuakes        <- nrow(quakes)
+
+  relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
+                               quakes[[numQuakes,  'ending']]), f=0.02)
+  threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+  plotMask   <- threeQuakeTimeMask
+  plotIndices<- which(plotMask)
+  plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
+
+  write(relaxedTime[[1]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'min', 'threequakes',
+                                    basedir), 'tex')))
+  write(relaxedTime[[2]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'max', 'threequakes',
+                                    basedir), 'tex')))
+
+  timeWindow = realTime[plotIndices]
+
+  relativeTau          <- h5file['relativeTimeIncrement'][]
+  fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
+  multiGridIterations  <- h5file['/iterations/multiGrid/final'][]
+  write.csv(data.frame(time = timeWindow,
+                       timeIncrement = finalTime * relativeTau[plotIndices],
+                       fixedPointIterations = fixedPointIterations[plotIndices],
+                       multiGridIterations = multiGridIterations[plotIndices]),
+            file.path(directories[['output']],
+                      paste.(pasteColon('2d-performance', basedir), 'csv')),
+            row.names = FALSE, quote = FALSE)
+  h5::h5close(h5file)
+}

data/tools/2d-velocity-contours.R

+source('tools/support/findQuakes.R')
+source('tools/support/writeContours.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+
+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (basedir in c("rfpitol=100e-7")) {
+  dir <- file.path(directories[['simulation']],
+                   '2d-lab-fpi-tolerance', basedir)
+  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+  relativeTime <- h5file['relativeTime'][]
+  realTime     <- finalTime * relativeTime
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]
+  basalTrenchDistance <- basalCoordinates[,1]
+  perm                <- order(basalTrenchDistance)
+  sortedBasalTrenchDistance <- basalTrenchDistance[perm]
+
+  {
+    ## We are interested in an enlarged time range around actual events here,
+    ## (and no other quantities!) hence we pass a very low velocity here.
+    quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
+                         indices = 1:dim(velocityProxy)[1], 1)
+    quakes$beginning <- realTime[quakes$beginningIndex]
+    quakes$ending    <- realTime[quakes$endingIndex]
+    quakes$duration  <- quakes$ending - quakes$beginning
+    numQuakes        <- nrow(quakes)
+
+    relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
+                                 quakes[[numQuakes,  'ending']]), f=0.02)
+    plotMask   <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+    plotIndices<- which(plotMask)
+
+    write(relaxedTime[[1]],
+          file.path(directories[['output']],
+                    paste.(pasteColon('timeframe', 'min', 'threequakes',
+                                      basedir), 'tex')))
+    write(relaxedTime[[2]],
+          file.path(directories[['output']],
+                    paste.(pasteColon('timeframe', 'max', 'threequakes',
+                                      basedir), 'tex')))
+
+    printlevels <- c('1000','100','10','1')
+    levels      <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
+    ret <- contourLines(realTime[plotIndices],
+                        sortedBasalTrenchDistance,
+                        abs(velocityProxy[plotIndices,perm,1][,,1]),
+                        levels = levels)
+
+    for (i in seq(printlevels))
+      writeContours(ret, levels[[i]],
+                    file.path(directories[['output']],
+                              paste.(pasteColon('2d-velocity-contours',
+                                                'threequakes', basedir, 'level',
+                                                printlevels[[i]]), 'tex')))
+  }
+  {
+    ## We are interested in an enlarged time range around actual events here,
+    ## (and no other quantities!) hence we pass a rather low velocity here.
+    quakes <- findQuakes(300e-6 + convergenceVelocity, velocityProxy,
+                         indices = 1:dim(velocityProxy)[1], 1)
+    quakes$beginning <- realTime[quakes$beginningIndex]
+    quakes$ending    <- realTime[quakes$endingIndex]
+    quakes$duration  <- quakes$ending - quakes$beginning
+    numQuakes        <- nrow(quakes)
+    quake            <- quakes[numQuakes,]
+    relaxedTime      <-
+      c(quake[['beginning']] - 0.9*(quake[['ending']] - quake[['beginning']]),
+        quake[['ending']]    + 0.1*(quake[['ending']] - quake[['beginning']]))
+    plotMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+    plotIndices<- which(plotMask)
+
+    printlevels <- c('3000','1000','300','100','30','10','3','1')
+    levels      <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
+    ret <- contourLines(realTime[plotIndices],
+                        sortedBasalTrenchDistance,
+                        abs(velocityProxy[plotIndices,perm,1][,,1]),
+                        levels = levels)
+
+    for (i in seq(printlevels))
+      writeContours(ret, levels[[i]],
+                    file.path(directories[['output']],
+                              paste.(pasteColon('2d-velocity-contours',
+                                                'zoom', basedir, 'level',
+                                                printlevels[[i]]), 'tex')))
+  }
+  h5::h5close(h5file)
+}

data/tools/3d-performance.R

+source('tools/support/findQuakes.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+
+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (basedir in c("rtol=1e-5_diam=1e-2")) {
+  dir          <- file.path(directories[['simulation']],
+                            '3d-lab', basedir)
+  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+  relativeTime <- h5file['relativeTime'][]
+  realTime     <- finalTime * relativeTime
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  ## We are interested in an enlarged time range around actual events here,
+  ## (and no other quantities!) hence we pass a very low velocity here.
+  quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
+                       ## Note: We only look at the last 2000 timesteps here because
+                       ## we're only interested in the last few events. This
+                       ## dramatically reduces RAM usage and runtime.
+                       indices = seq(dim(velocityProxy)[1]-2000,
+                                     dim(velocityProxy)[1]), c(1,3))
+  quakes$beginning <- realTime[quakes$beginningIndex]
+  quakes$ending    <- realTime[quakes$endingIndex]
+  quakes$duration  <- quakes$ending - quakes$beginning
+  numQuakes        <- nrow(quakes)
+
+  relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
+                               quakes[[numQuakes,  'ending']]), f=0.02)
+  threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
+  plotMask   <- threeQuakeTimeMask
+  plotIndices<- which(plotMask)
+  plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)
+
+  write(relaxedTime[[1]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'min', 'threequakes',
+                                    basedir), 'tex')))
+  write(relaxedTime[[2]],
+        file.path(directories[['output']],
+                  paste.(pasteColon('timeframe', 'max', 'threequakes',
+                                    basedir), 'tex')))
+
+  timeWindow = realTime[plotIndices]
+
+  relativeTau          <- h5file['relativeTimeIncrement'][]
+  fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
+  multiGridIterations  <- h5file['/iterations/multiGrid/final'][]
+  write.csv(data.frame(time = timeWindow,
+                       timeIncrement = finalTime * relativeTau[plotIndices],
+                       fixedPointIterations = fixedPointIterations[plotIndices],
+                       multiGridIterations = multiGridIterations[plotIndices]),
+            file.path(directories[['output']],
+                      paste.(pasteColon('3d-performance', basedir), 'csv')),
+            row.names = FALSE, quote = FALSE)
+  h5::h5close(h5file)
+}

data/tools/3d-velocity-contours.R

+source('tools/support/findQuakes.R')
+source('tools/support/writeContours.R')
+
+finalTime           <- 1000 # s
+convergenceVelocity <- 5e-5 # m/s
+printlevels <- c('1','2','3','5',
+                 '10','20','30','50',
+                 '100','200','300','500',
+                 '1000')
+criticalVelocities  <- 1e-6*as.numeric(printlevels) + convergenceVelocity
+
+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+directories <- ini::read.ini('config.ini')$directories
+dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
+for (basedir in c("rtol=1e-5_diam=1e-2")) {
+  dir          <- file.path(directories[['simulation']],
+                            '3d-lab', basedir)
+  h5file       <- h5::h5file(file.path(dir, 'output.h5'), 'r')
+  relativeTime <- h5file['relativeTime'][]
+  realTime     <- finalTime * relativeTime
+
+  velocityProxy<- h5file['/frictionalBoundary/velocity']
+
+  ## We are interested in an enlarged time range around actual events here,
+  ## (and no other quantities!) hence we pass a rather low velocity here.
+  quakes <- findQuakes(200e-6 + convergenceVelocity, velocityProxy,
+                       ## Note: We only look at the last 1000 timesteps here because
+                       ## we're only interested in the last few events. This
+                       ## dramatically reduces RAM usage and runtime.
+                       indices = seq(dim(velocityProxy)[1]-1000,
+                                     dim(velocityProxy)[1]), c(1,3))
+  quake <- quakes[nrow(quakes)-1,] # Q: why did we not need the -1 in julia?
+  
+  weakPatchGridVelocityProxy <- h5file["/weakPatchGrid/velocity"]
+
+  stepSize <- 30 ## note: should/might differ by time/spatial resolution
+  ts <- seq(quake$beginningIndex, quake$endingIndex, by=stepSize)
+  dd <- data.frame(timeSteps = ts,
+                   times = realTime[ts],
+                   timeOffsets = realTime[ts] - realTime[quake$beginningIndex])
+
+  fname = pasteColon('3d-velocity-contours', basedir)
+  write.csv(dd, row.names = FALSE, quote = FALSE,
+            file = file.path(directories[['output']],
+                             paste.(pasteColon(fname, 'times'), 'csv')))
+
+  weakPatchGridXCoordinates <- h5file["/weakPatchGrid/xCoordinates"][]
+  weakPatchGridZCoordinates <- h5file["/weakPatchGrid/zCoordinates"][]
+  
+  for (t in ts) {
+    m <- weakPatchGridVelocityProxy[t,,,]
+    s <- sqrt(m[,,,1]^2 + m[,,,3]^2)
+    ret <- contourLines(weakPatchGridXCoordinates, weakPatchGridZCoordinates, s,
+                        level=criticalVelocities)
+
+    for (i in seq(printlevels))
+      writeContours(ret, criticalVelocities[[i]],
+                    file.path(directories[['output']],
+                              paste.(pasteColon(fname,
+                                                'step', t,
+                                                'level', printlevels[[i]]),
+                                     'tex')))
+  }
+  h5::h5close(h5file)
+}

data/tools/comparison:lab-sim.R

+paste.    <- function(...) paste(..., sep='.')
+pasteColon<- function(...) paste(..., sep=':')
+
+basedir <- 'rfpitol=100e-7'
+
+directories <- ini::read.ini('config.ini')$directories
+labdata <- within(
+  read.table(unz(file.path(directories[['experiment']],
+                           'B_Scale-model-earthquake-data.zip'),
+                 'scaleEQs_12-31_catalogue.txt'),
+             sep='\t', header=FALSE, comment.char='%',
+             col.names=c('frame',
+                         'meanSlipNature', 'meanSlipLab',
+                         'peakSlipNature', 'peakSlipLab',
+                         'ruptureWidthNature', 'ruptureWidthLab',
+                         'ignored', 'ignored', 'ignored', 'ignored', 'ignored')),
+  {
+    recurrenceLab <- c(NA, diff(frame))/10 # (10Hz camera: 10 frames ~ 1s)
+    meanSlipLab   <- meanSlipLab / 1e6     # micro m -> m
+    peakSlipLab   <- peakSlipLab / 1e6})   # micro m -> m
+## NB: We skip the first two quakes here, for compatibility with my old code
+##     Maybe skipping the first was intentional so that each quake would have
+##     a recurrence time. But skipping the second was certainly an accident
+labdata <- labdata[3:nrow(labdata),]
+
+simdata <- read.csv(file.path(directories[['output']],
+                              paste.(pasteColon('events', basedir), 'csv')))
+
+report <- function(lab, sim, quantity) {
+  write.table(lab, file.path(directories[['output']],
+                             paste.(pasteColon('boxplot-data', 'lab',
+                                               quantity), 'tex')),
+              row.names=FALSE, col.names=FALSE)
+  write.table(sim, file.path(directories[['output']],
+                             paste.(pasteColon('boxplot-data', 'simulation',
+                                               basedir, quantity), 'tex')),
+              row.names=FALSE, col.names=FALSE)
+}
+report(labdata$recurrenceLab, diff(simdata$beginning), 'recurrence')
+report(labdata$ruptureWidthLab, simdata$ruptureWidth, 'ruptureWidth')
+# meters -> millimeters
+report(labdata$peakSlipLab * 1000, simdata$peakSlip*1000, 'peakSlip')

data/tools/generate-2d.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# contours
+rr 2d-velocity-contours.R
+
+# dip
+rr 2d-dip-contours.R
+
+# iterations / adaptivity
+rr 2d-performance.R
+
+# box plot (one half)
+rr 2d-event-writer.R
+
+# fpi data
+rr 2d-fpi-tolerance.R
+
+date

data/tools/generate-3d.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# performance
+rr 3d-performance.R
+
+# velocity contours
+rr 3d-velocity-contours.R
+
+date

data/tools/generate-others.bash

+#!/usr/bin/env bash
+
+set -e
+
+rr() {
+    echo "$(date)" Running: Rscript $@
+    Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
+}
+
+# box plot (one half)
+rr comparison:lab-sim.R
+
+date

data/tools/support/findQuakes.R

+Rcpp::sourceCpp('tools/support/maxVelocity.cpp')
+Rcpp::sourceCpp('tools/support/negativeStarts.cpp')
+Rcpp::sourceCpp('tools/support/positiveStarts.cpp')
+
+findQuakes <- function (criticalVelocity, velocity, indices, dimensions) {
+  maximumVelocities<- maxVelocity(velocity[indices,,dimensions])
+  slippingTimes  <- maximumVelocities > criticalVelocity
+
+  ns <- negativeStarts(slippingTimes) + indices[[1]] - 1
+  ps <- positiveStarts(slippingTimes) + indices[[1]] - 1
+
+  ## NOTE: we drop incomplete quakes here
+  mlen <- min(length(ns), length(ps))
+  ns <- ns[1:mlen]
+  ps <- ps[1:mlen]
+
+  return(data.frame(endingIndex = ns, beginningIndex = ps))
+}

data/tools/support/maxVelocity.cpp

+#include <Rcpp.h>
+
+// For std::hypot(x,y)
+// [[Rcpp::plugins(cpp11)]]
+
+// Layout of vectors is such that loops over inner indices are traversed first
+// [[Rcpp::export]]
+Rcpp::NumericVector maxVelocity(Rcpp::NumericVector const &x) {
+  Rcpp::IntegerVector size = x.attr("dim");
+  Rcpp::NumericVector ret(size[0]);
+  switch (size[2]) {
+  case 1:
+    for (size_t ts = 0; ts < size[0]; ++ts)
+      for (size_t coord = 0; coord < size[1]; ++coord)
+        ret[ts] = std::max(ret[ts],
+                           std::abs(x[0 * size[1] * size[0] + coord * size[0] + ts]));
+    break;
+  case 2:
+    for (size_t ts = 0; ts < size[0]; ++ts)
+      for (size_t coord = 0; coord < size[1]; ++coord)
+        ret[ts] = std::max(ret[ts],
+                           std::hypot(x[0 * size[1] * size[0] + coord * size[0] + ts],
+                                      x[1 * size[1] * size[0] + coord * size[0] + ts]));
+    break;
+  default:
+    throw std::range_error("Inadmissible value");
+  }
+  return ret;
+}

data/tools/support/negativeStarts.cpp

+#include <Rcpp.h>
+
+// [[Rcpp::export]]
+Rcpp::NumericVector negativeStarts(Rcpp::NumericVector const &x) {
+  std::vector<size_t> starts(0);
+  bool prev = false;
+  for (size_t i = 0; i < x.size(); ++i) {
+    if (prev && !x[i])
+      starts.push_back(i+1); // R indices start at 1
+    prev = x[i];
+  }
+  return Rcpp::NumericVector(starts.begin(), starts.end());
+}

data/tools/support/positiveStarts.cpp

+#include <Rcpp.h>
+
+// [[Rcpp::export]]
+Rcpp::NumericVector positiveStarts(Rcpp::NumericVector const &x) {
+  std::vector<size_t> starts(0);
+  bool prev = false;
+  for (size_t i = 0; i < x.size(); ++i) {
+    if (!prev && x[i])
+      starts.push_back(i+1); // R indices start at 1
+    prev = x[i];
+  }
+  return Rcpp::NumericVector(starts.begin(), starts.end());
+}
+

data/tools/support/trapezoidal.cpp

+#include <Rcpp.h>
+
+// [[Rcpp::export]]
+double trapezoidal(Rcpp::NumericVector const &x, Rcpp::NumericVector const &y) {
+  double ret = 0;
+  for (size_t i = 1; i < x.size(); ++i)
+    ret += (x[i] - x[i - 1]) * (y[i] + y[i - 1]) / 2;
+  return ret;
+}

data/tools/support/writeContours.R

+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)
+  }
+}

debugging.m

+
+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