diff --git a/receiver/bin/viewrec b/receiver/bin/viewrec new file mode 100755 index 0000000..ec0557f --- /dev/null +++ b/receiver/bin/viewrec @@ -0,0 +1,3 @@ +#!/bin/bash +DIR=$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd ) +python3 $DIR/../src/viewrec.py diff --git a/receiver/src/Filters.py b/receiver/src/Filters.py new file mode 100644 index 0000000..804cdf2 --- /dev/null +++ b/receiver/src/Filters.py @@ -0,0 +1,244 @@ +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2016, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +from PyQt5.QtCore import * +from PyQt5.QtWidgets import * +import math +import scipy.signal +import numpy +import abc + +class Filter(QGroupBox): + activeItemChanged = pyqtSignal(name='filterChanged') + + def __init__(self, title, parent = None): + super(Filter, self).__init__(title, parent) + + self.setCheckable(True) + self.setChecked(False) + + self.toggled.connect(self.filterChanged) + + @abc.abstractmethod + def apply(self, wf): + """ Filter waveforms. """ + return + +class Lowpass(Filter): + def __init__(self, parent = None): + super(Lowpass, self).__init__('Lowpass', parent) + + lowpassOrderLabel = QLabel('Order', self) + self.lowpassOrder = QSpinBox(self) + self.lowpassOrder.setValue(4) + self.lowpassOrder.valueChanged.connect(self.filterChanged) + cutoffLabel = QLabel('Cutoff (Hz)', self) + self.cutoff = QDoubleSpinBox(self) + self.cutoff.setValue(0.5) + self.cutoff.setSingleStep(0.1) + self.cutoff.valueChanged.connect(self.filterChanged) + + filterLayout = QFormLayout(self) + filterLayout.addRow(lowpassOrderLabel, self.lowpassOrder) + filterLayout.addRow(cutoffLabel, self.cutoff) + + def apply(self, wf): + Fs = 1.0 / (wf.time[1] - wf.time[0]) + cutoff = min(0.5 * Fs, max(1. / Fs, self.cutoff.value())) + b, a = scipy.signal.butter(self.lowpassOrder.value(), cutoff, 'low', fs=Fs) + for name in wf.waveforms.keys(): + wf.waveforms[name] = scipy.signal.filtfilt(b, a, wf.waveforms[name]) + +class Deconvolve(Filter): + def __init__(self, parent = None): + super(Deconvolve, self).__init__('Deconvolve', parent) + + inLabel = QLabel('Input', self) + inFun = QLabel('t/T^2 * exp(-t/T)', self) + TLabel = QLabel('T', self) + self.T = QDoubleSpinBox(self) + self.T.setValue(0.1) + self.T.setMaximum(float('inf')) + self.T.valueChanged.connect(self.filterChanged) + + outLabel = QLabel('Output', self) + outFun = QLabel('Gauss', self) + sigmaLabel = QLabel('Sigma', self) + self.sigma = QDoubleSpinBox(self) + self.sigma.setValue(0.05) + self.sigma.setMaximum(float('inf')) + self.sigma.valueChanged.connect(self.filterChanged) + + filterLayout = QFormLayout(self) + filterLayout.addRow(inLabel, inFun) + filterLayout.addRow(TLabel, self.T) + filterLayout.addRow(outLabel, outFun) + filterLayout.addRow(sigmaLabel, self.sigma) + + def deconv(self, waveform, dt): + y = numpy.zeros(waveform.size) + T = self.T.value() + s = self.sigma.value() + F = lambda t : numpy.exp(-numpy.square(t-4*s) / (2*s**2)) / (math.sqrt(2*math.pi) * s) + G = lambda t : 1 - 2*T/s**2 * (t-4*s) - T**2/s**2 * (1 - numpy.square(t-4*s) / s**2) + eps = 1e-16 + th = s * (4 + math.sqrt(-2 * math.log(eps))) # exp(-(t-4*s)^2 / (2*s**2)) is smaller than eps if t >= th + + # convolves waveform with F*T. The integral int f(t) dt is approximated as simple and stupid as sum_i f(t_i) * dt + for j in range(y.size): + i = numpy.arange(max(0, int(j-th/dt)), j) + t = (j - i) * dt + f = numpy.multiply(F(t), G(t)) + y[j] = numpy.sum(dt * numpy.multiply(waveform[i], f)) + + return y + + + def apply(self, wf): + dt = wf.time[1] - wf.time[0] + keys = [f'v{i+1}' for i in range(3)] + for k in keys: + if k in wf.waveforms: + wf.waveforms[k] = self.deconv(wf.waveforms[k], dt) + +class Rotate(Filter): + def __init__(self, parent = None): + super(Rotate, self).__init__('Rotate', parent) + + coordsysLabel = QLabel('Coordinates', self) + self.coordsys = QComboBox(self) + self.coordsys.addItem('ned') + self.coordsys.addItem('seu') + self.coordsys.currentIndexChanged.connect(self.filterChanged) + + epicenterXLabel = QLabel('Epicenter X', self) + self.epicenterX = QDoubleSpinBox(self) + self.epicenterX.setValue(0.0) + self.epicenterX.setMaximum(float('inf')) + self.epicenterX.valueChanged.connect(self.filterChanged) + epicenterYLabel = QLabel('Epicenter Y', self) + self.epicenterY = QDoubleSpinBox(self) + self.epicenterY.setValue(0.0) + self.epicenterY.setMaximum(float('inf')) + self.epicenterY.valueChanged.connect(self.filterChanged) + + filterLayout = QFormLayout(self) + filterLayout.addRow(coordsysLabel, self.coordsys) + filterLayout.addRow(epicenterXLabel, self.epicenterX) + filterLayout.addRow(epicenterYLabel, self.epicenterY) + + def apply(self, wf): + if 'v1' in wf.waveforms and 'v2' in wf.waveforms and 'v3' in wf.waveforms: + epicenter = numpy.array([self.epicenterX.value(), self.epicenterY.value(), 0.0]) + radial = wf.coordinates - epicenter + phi = math.acos(radial[0] / numpy.linalg.norm(radial)) + + u = wf.waveforms.pop('v1') + v = wf.waveforms.pop('v2') + w = wf.waveforms.pop('v3') + + if self.coordsys.currentText() == 'seu': + u = -u + w = -w + + wf.waveforms['radial'] = math.cos(phi) * u + math.sin(phi) * v + wf.waveforms['transverse'] = -math.sin(phi) * u + math.cos(phi) * v + wf.waveforms['vertical'] = w + for new_comp in ['radial', 'transverse', 'vertical']: + wf.show[new_comp] = True + +class Pick(Filter): + def __init__(self, parent = None): + super(Pick, self).__init__('Pick components', parent) + self.cb_widget_list = [] + self.layout = QGridLayout() + self.setLayout(self.layout) + + def create_checkbox(self, name): + widget = QCheckBox(name) + widget.stateChanged.connect(self.filterChanged) + return widget + + def add_widgts_to_layout(self): + cols = 3 + nwidget_over_3 = len(self.cb_widget_list) // cols + rows = nwidget_over_3 + 1 + for k, widget in enumerate(self.cb_widget_list): + row = k // cols + col = k - row * cols + self.layout.addWidget(widget, row, col) + + def create_checkboxes(self, wf): + assert self.cb_widget_list == [] + + for name in wf.waveforms.keys(): + self.cb_widget_list.append(self.create_checkbox(name)) + + self.add_widgts_to_layout() + + def update_checkboxes(self, wf): + existing_widgets = [cb.text() for cb in self.cb_widget_list] + required_widgets = list(wf.waveforms.keys()) + + # remove not needed checkboxes + # can't remove from list, while iterating over list + widgets_to_remove = [] + for w in self.cb_widget_list: + if not w.text() in required_widgets: + widgets_to_remove.append(w) + for w in widgets_to_remove: + w.deleteLater() + self.cb_widget_list.remove(w) + + # add new checkboxes + for name in required_widgets: + if not name in existing_widgets: + self.cb_widget_list.append(self.create_checkbox(name)) + + self.add_widgts_to_layout() + + def apply(self, wf): + if not self.cb_widget_list: + self.create_checkboxes(wf) + else: + self.update_checkboxes(wf) + + for widget in self.cb_widget_list: + var_name = widget.text() + wf.show[var_name] = widget.isChecked() diff --git a/receiver/src/Navigation.py b/receiver/src/Navigation.py new file mode 100644 index 0000000..e4941ad --- /dev/null +++ b/receiver/src/Navigation.py @@ -0,0 +1,142 @@ +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger) +# +# @section LICENSE +# Copyright (c) 2015-2016, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +from PyQt5.QtCore import * +from PyQt5.QtGui import * +from PyQt5.QtWidgets import * + +import os +import copy + +import Tecplot +import Waveform + +class Navigation(QWidget): + activeItemChanged = pyqtSignal(name='activeItemChanged') + folderChanged = pyqtSignal(str, str, name='folderChanged') + close = pyqtSignal(QWidget, name='close') + + def __init__(self, noclose = False, parent = None): + super(Navigation, self).__init__(parent) + + self.currentFolder = '' + + openIcon = QIcon.fromTheme('folder-open') + openButton = QPushButton(openIcon, '', self) + openButton.clicked.connect(self.selectFolder) + refreshIcon = QIcon.fromTheme('view-refresh') + refreshButton = QPushButton(refreshIcon, '', self) + refreshButton.clicked.connect(self.refreshFolder) + if not noclose: + closeIcon = QIcon.fromTheme('window-close') + closeButton = QPushButton(closeIcon, '', self) + closeButton.clicked.connect(self.emitClose) + + self.receiverList = QListView(self) + self.model = QStandardItemModel() + self.receiverList.setModel(self.model) + self.receiverList.clicked.connect(self.activeItemChanged) + + buttonLayout = QHBoxLayout() + buttonLayout.addWidget(openButton) + buttonLayout.addWidget(refreshButton) + if not noclose: + buttonLayout.addWidget(closeButton) + buttonLayout.addStretch() + + layout = QVBoxLayout(self) + layout.addLayout(buttonLayout) + self.receiverList.setMaximumWidth(120) + layout.addWidget(self.receiverList) + + def selectFolder(self): + folder = QFileDialog.getExistingDirectory(self, 'Open directory', self.currentFolder, QFileDialog.ShowDirsOnly) + if not folder=='': + folder = str(folder) + self.readFolder(folder) + self.folderChanged.emit(self.currentFolder, folder) + self.currentFolder = folder + + def readFolder(self, folder): + if len(folder) != 0: + currentIndex = self.receiverList.currentIndex() + self.model.clear() + + files = [f for f in os.listdir(folder) if os.path.isfile(os.path.join(folder,f))] + files.sort() + numValidFiles = 0 + for f in files: + if f.endswith('dat'): + wf = Tecplot.read(os.path.join(folder,f)) + if wf: + item = QStandardItem(f) + item.setData(wf) + self.model.appendRow(item) + numValidFiles = numValidFiles + 1 + + if currentIndex.row() >= 0 and numValidFiles > currentIndex.row(): + newIndex = self.model.index(currentIndex.row(), currentIndex.column()) + self.receiverList.setCurrentIndex(newIndex) + self.activeItemChanged.emit() + + def getActiveWaveforms(self): + waveforms = [] + for index in self.receiverList.selectedIndexes(): + wf = self.model.itemFromIndex(index).data() + waveforms.append( copy.deepcopy(wf) ) + return waveforms + + def refreshFolder(self): + self.readFolder(self.currentFolder) + if not self.receiverList.selectionModel().hasSelection(): + self.activeItemChanged.emit() + + def numberOfRows(self): + return self.model.rowCount() + + def selectWaveformAt(self, row): + self.receiverList.selectionModel().select(self.model.index(row, 0), QItemSelectionModel.ClearAndSelect) + + def emitClose(self): + self.folderChanged.emit(self.currentFolder, '') + self.close.emit(self) + + + diff --git a/receiver/src/Tecplot.py b/receiver/src/Tecplot.py new file mode 100644 index 0000000..81c400f --- /dev/null +++ b/receiver/src/Tecplot.py @@ -0,0 +1,83 @@ +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2015, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +import Waveform +import re + +def read(fileName): + data = [] + coordinates = [float('nan')] * 3 + coordComment = re.compile('#\s*x(\d)\s+([0-9\.eE\+\-]+)') + variables = '' + p_0 = 0. + t_s = 0. + t_d = 0. + f = open(fileName) + for row in f: + row = row.strip() + if not row[0].isalpha() and not row[0] == '#': + values = row.split() + data.append([float(numeric_string) for numeric_string in values]) + elif row[0] == '#': + match = coordComment.match(row) + if match: + coordinates[int(match.group(1))-1] = float(match.group(2)) + elif row[2:5] == 'P_0': + p_0 = float(row[5:]) + elif row[2:5] == 'T_s': + t_s = float(row[5:]) + elif row[2:5] == 'T_d': + t_d = float(row[5:]) + elif row.startswith('VARIABLES'): + row = row.split('=') + variables = row[1].replace('"', '').split(',') + variables = [v.strip() for v in variables] + f.close() + if len(data) == 0: + return None + if 'P_n' in variables: + indexP_n = variables.index('P_n') + data = [d[0:indexP_n] + [d[indexP_n] + p_0] + d[indexP_n+1:] for d in data] + if 'T_s' in variables: + indexT_s = variables.index('T_s') + data = [d[0:indexT_s] + [d[indexT_s] + t_s] + d[indexT_s+1:] for d in data] + if 'T_d' in variables: + indexT_d = variables.index('T_d') + data = [d[0:indexT_d] + [d[indexT_d] + t_d] + d[indexT_d+1:] for d in data] + return Waveform.Waveform(variables, data, coordinates) diff --git a/receiver/src/View.py b/receiver/src/View.py new file mode 100644 index 0000000..89f46c4 --- /dev/null +++ b/receiver/src/View.py @@ -0,0 +1,258 @@ +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger) +# +# @section LICENSE +# Copyright (c) 2015-2016, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +from PyQt5.QtGui import * +from PyQt5.QtWidgets import * +from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas +try: + from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT as NavigationToolbar +except ImportError: + from matplotlib.backends.backend_qt5agg import NavigationToolbar2QTAgg as NavigationToolbar +import matplotlib.pyplot as plt + +import Navigation +import Filters +import Watchdog +import re +import math +import numpy +import os.path +import scipy.fftpack +import scipy.interpolate +import scipy.integrate + +class View(QWidget): + + def __init__(self, parent = None): + super(View, self).__init__(parent) + + self.__watchdog = Watchdog.Watchdog() + self.__watchdog.fileChanged.connect(self.refreshAll) + + self.figure = plt.figure() + self.canvas = FigureCanvas(self.figure) + self.canvas.setSizePolicy(QSizePolicy.Expanding, QSizePolicy.Expanding) + toolbar = NavigationToolbar(self.canvas, self) + + self.navigationLayout = QHBoxLayout() + + layout = QHBoxLayout(self) + self.navigations = [] + self.addNavigation(True) + + self.filters = [ Filters.Lowpass(), Filters.Deconvolve(), Filters.Rotate(), Filters.Pick() ] + filterLayout = QVBoxLayout() + for f in self.filters: + filterLayout.addWidget(f) + f.filterChanged.connect(self.plot) + filterLayout.addStretch() + + addIcon = QIcon.fromTheme('list-add') + addNaviButton = QPushButton(addIcon, 'Add navigation', self) + addNaviButton.clicked.connect(self.addNavigation) + + self.maxFreq = QDoubleSpinBox(self) + self.maxFreq.setValue(10.0) + self.maxFreq.setVisible(False) + self.maxFreq.valueChanged.connect(self.plot) + spectrumIcon = QIcon.fromTheme('network-wireless') + self.spectrum = QPushButton(spectrumIcon, 'Spectrum', self) + self.spectrum.setCheckable(True) + self.spectrum.clicked.connect(self.plot) + self.spectrum.toggled.connect(self.maxFreq.setVisible) + self.diff = QPushButton('Diff', self) + self.diff.setCheckable(True) + self.diff.clicked.connect(self.plot) + self.diff.clicked.connect(self.spectrum.setHidden) + self.spectrum.toggled.connect(self.diff.setHidden) + + autoRefresh = QPushButton(QIcon.fromTheme('view-refresh'), 'Auto', self) + autoRefresh.setCheckable(True) + autoRefresh.clicked.connect(self.__watchdog.toggle) + + saveAll = QPushButton(QIcon.fromTheme('document-save'), '', self) + saveAll.clicked.connect(self.savePlots) + + self.normalize = QPushButton('Normalize', self) + self.normalize.setCheckable(True) + self.normalize.clicked.connect(self.plot) + + self.differentiate = QPushButton('Differentiate', self) + self.differentiate.setCheckable(True) + self.differentiate.clicked.connect(self.plot) + + toolLayout = QHBoxLayout() + toolLayout.addWidget(addNaviButton) + toolLayout.addWidget(self.diff) + toolLayout.addWidget(self.spectrum) + toolLayout.addWidget(self.maxFreq) + toolLayout.addWidget(autoRefresh) + toolLayout.addWidget(saveAll) + toolLayout.addWidget(self.normalize) + toolLayout.addWidget(self.differentiate) + toolLayout.addWidget(toolbar) + plotLayout = QVBoxLayout() + plotLayout.addLayout(toolLayout) + plotLayout.addWidget(self.canvas) + layout.addLayout(self.navigationLayout) + layout.addLayout(plotLayout) + layout.addLayout(filterLayout) + + def addNavigation(self, noclose = False): + navigation = Navigation.Navigation(noclose) + navigation.activeItemChanged.connect(self.plot) + navigation.folderChanged.connect(self.navigationFolderChanged) + navigation.close.connect(self.closeNavigation) + navigation.setSizePolicy(QSizePolicy.Maximum, QSizePolicy.Minimum) + self.navigationLayout.addWidget(navigation) + self.navigations.append(navigation) + + def navigationFolderChanged(self, oldFolder, newFolder): + self.__watchdog.removeFolder(oldFolder) + self.__watchdog.addFolder(newFolder) + + def closeNavigation(self, widget): + self.navigations.remove(widget) + self.navigationLayout.removeWidget(widget) + widget.deleteLater() + self.plot() + + def refreshAll(self): + for navigation in self.navigations: + navigation.refreshFolder() + + def plot(self): + wfc = [wf for nav in self.navigations for wf in nav.getActiveWaveforms()] + for filt in self.filters: + if filt.isChecked(): + for wf in wfc: + filt.apply(wf) + + if self.normalize.isChecked(): + #normalize traces + for nWf, wf in enumerate(wfc): + wfc[nWf].normalize() + + if self.differentiate.isChecked(): + #differentiate traces + for nWf, wf in enumerate(wfc): + wfc[nWf].differentiate() + + if self.diff.isChecked() and len(wfc) > 0: + wf0 = wfc.pop() + for nWf, wf in enumerate(wfc): + wfc[nWf].subtract(wf0) + + names = set([name for wf in wfc for name in wf.waveforms.keys() if wf.show[name]]) + numPlots = len(names) + + self.figure.clear() + if numPlots > 0: + names = list(names) + names.sort() + + numRows = math.ceil(math.sqrt(numPlots)); + numCols = math.ceil(numPlots / numRows) + subplots = dict() + for i in range(numPlots): + subplots[ names[i] ] = self.figure.add_subplot(numRows, numCols, i+1) + + wf_ref = wfc[0] + for nWf, wf in enumerate(wfc): + for name, waveform in wf.waveforms.items(): + if not wf.show[name]: + continue + p = subplots[name] + if self.spectrum.isChecked(): + n = len(waveform) + dt = wf.time[1]-wf.time[0] # assume equally spaced samples + f = scipy.fftpack.fftfreq(n, dt) + W = dt * scipy.fftpack.fft(waveform) + maxFreqIndices = numpy.argwhere(f > self.maxFreq.value()) + L = maxFreqIndices[0,0] if numpy.size(maxFreqIndices) > 0 else n/2 + p.loglog(f[1:L], numpy.absolute(W[1:L]), label=str(nWf)) + p.set_xlabel('f [Hz]') + elif self.diff.isChecked(): + p.plot(wf.time, waveform, label='{}-0'.format(nWf+1)) + p.set_xlabel('t (s)') + else: + p.plot(wf.time, waveform, label=str(nWf)) + p.set_xlabel('t (s)') + p.set_ylabel(name) + #print L2 difference + if nWf > 0 and not self.diff.isChecked() and name in wf_ref.waveforms: + t_min = max(wf.time.min(), wf_ref.time.min()) + t_max = min(wf.time.max(), wf_ref.time.max()) + truncate = lambda a: a[(a >= t_min) & (a <= t_max)] + time_union = numpy.union1d(truncate(wf.time), truncate(wf_ref.time)) + wf_interp = scipy.interpolate.interp1d(wf.time, waveform) + ref_interp = scipy.interpolate.interp1d(wf_ref.time, wf_ref.waveforms[name]) + wf_union = wf_interp(time_union) + ref_union = ref_interp(time_union) + diff = numpy.sqrt(scipy.integrate.trapz((wf_union - ref_union)**2, x=time_union)) + p.text(0.05, 1-0.1*nWf, f"L2 diff={diff:.2e}", transform = p.transAxes) + + self.figure.tight_layout() + + for i in range(len(names)): + if subplots[ names[i] ]: + subplots[ names[i] ].legend(prop={'size':8}, frameon=False) + self.canvas.draw() + + def savePlots(self): + filetypes = self.canvas.get_supported_filetypes_grouped() + defaultFiletype = self.canvas.get_default_filetype() + filters = [] + selectedFilter = '' + for name, extensions in sorted(filetypes.items()): + filtr = '{0} ({1})'.format(name, ' '.join(['*.{0}'.format(ext) for ext in extensions])) + if defaultFiletype in extensions: + selectedFilter = filtr + filters.append(filtr) + fileName, filtr = QFileDialog.getSaveFileName(self, 'Choose a save location.', '', ';;'.join(filters), selectedFilter) + fileName = os.path.splitext(str(fileName))[0] + extension = re.search(r'\*(\.[a-zA-Z]+)', str(filtr)).group(1) + + maxRow = min([nav.numberOfRows() for nav in self.navigations]) + for row in range(maxRow): + for nav in self.navigations: + nav.selectWaveformAt(row) + self.plot() + self.canvas.print_figure('{0}{1:03}{2}'.format(fileName, row+1, extension)) diff --git a/receiver/src/Watchdog.py b/receiver/src/Watchdog.py new file mode 100644 index 0000000..7ee2a3d --- /dev/null +++ b/receiver/src/Watchdog.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger) +# +# @section LICENSE +# Copyright (c) 2016, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +from PyQt5.QtCore import QThread, pyqtSignal + +class Watchdog(QThread): + fileChanged = pyqtSignal(name='fileChanged') + + def __init__(self, parent = None): + QThread.__init__(self, parent) + + self.__enabled = False + self.__folders = [] + + def __del__(self): + self.__enabled = False + self.wait() + + def toggle(self): + self.__enabled = not self.__enabled + if self.__enabled: + self.start() + + def addFolder(self, folder): + if folder: + self.__folders.append(folder) + + def removeFolder(self, folder): + try: + self.__folders.remove(folder) + except ValueError: + # Can happen if folder was empty + pass + + def run(self): + while self.__enabled: + # TODO check files + self.fileChanged.emit() + self.sleep(3) diff --git a/receiver/src/Waveform.py b/receiver/src/Waveform.py new file mode 100644 index 0000000..9e9c35a --- /dev/null +++ b/receiver/src/Waveform.py @@ -0,0 +1,77 @@ +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2015, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +import numpy + +class Waveform: + def __init__(self, names, data, coordinates): + data = numpy.array(data) + + self.waveforms = dict() + self.norm = dict() + self.show = dict() + + for i in range(0, len(names)): + if names[i] == 'Time': + self.time = data[:,i] + else: + self.waveforms[ names[i] ] = data[:,i] + self.norm[ names[i] ] = numpy.max(numpy.abs(data[:,i])) + self.show[ names[i] ] = True + + self.coordinates = numpy.array(coordinates) + + def subtract(self, other): + newTime = numpy.union1d(self.time, other.time) + self.waveforms = { key: value for key,value in self.waveforms.items() if key in other.waveforms } + for name, wf in self.waveforms.items(): + wf0 = numpy.interp(newTime, other.time, other.waveforms[name]) + wf1 = numpy.interp(newTime, self.time, self.waveforms[name]) + self.waveforms[name] = wf1 - wf0 + self.time = newTime + + def normalize(self): + for name, wf in self.waveforms.items(): + if self.norm[name] > numpy.finfo(float).eps: + self.waveforms[name] = wf / self.norm[name] + + def differentiate(self): + dt = self.time[1] - self.time[0] + for name, wf in self.waveforms.items(): + self.waveforms[name] = numpy.gradient(wf, dt) diff --git a/receiver/src/__init__.py b/receiver/src/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/receiver/src/__init__.py @@ -0,0 +1 @@ + diff --git a/receiver/src/viewrec.py b/receiver/src/viewrec.py new file mode 100644 index 0000000..895a802 --- /dev/null +++ b/receiver/src/viewrec.py @@ -0,0 +1,52 @@ +#!/usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# +# @section LICENSE +# Copyright (c) 2015, SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# + +import sys +from PyQt5 import QtWidgets +import View + +app = QtWidgets.QApplication(sys.argv) + +w = View.View() +w.setWindowTitle('SeisSol receiver view') +w.resize(800, 600) +w.show() + +app.exec_() diff --git a/tools/CreateVtkCoastLineFromGmt.py b/tools/CreateVtkCoastLineFromGmt.py new file mode 100755 index 0000000..b29853c --- /dev/null +++ b/tools/CreateVtkCoastLineFromGmt.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +import argparse +import os +import numpy as np + +parser = argparse.ArgumentParser(description="create coastline vtk file from gmt") +parser.add_argument("--lon", nargs=2, metavar=(("lonmin"), ("lonmax")), help="lonmin: minimum longitude, lonmax: maximum longitude", type=float) +parser.add_argument("--lat", nargs=2, metavar=(("latmin"), ("latmax")), help="latmin: minimum latitude, lonmax: maximum latitude", type=float) +parser.add_argument("--proj", nargs=1, metavar=("projname"), help="project the data. projname: PROJ string describing the projection (ex epsg:32646 for UTM46N). Use geocent for geocentric cartesian") +parser.add_argument("--resolution", nargs=1, metavar=("resolution"), default=("i"), help="Coastline resolution, (f)ull, (h)igh, (i)ntermediate, (l)ow, and (c)rude") +parser.add_argument("--recenter", nargs=2, metavar=("x", "y"), default=([0, 0]), help="translate coordinate array (e.g. x_new = x_old - x)", type=float) +parser.add_argument("--z", nargs=1, metavar=("elevation"), default=([0]), help="z coordinate of coastline", type=float) +args = parser.parse_args() + + +# export cordinates from GMT +os.system("module load gmt") +os.system(f"gmt pscoast -R{args.lon[0]}/{args.lon[1]}/{args.lat[0]}/{args.lat[1]} -D{args.resolution[0]} -M -W > coastline.dat") + +# Read GMT file +xyz = [] +segments = [] +nvert = 0 +newPolyLine = True +with open("coastline.dat") as fid: + for line in fid: + if line.startswith("#"): + continue + if line.startswith(">"): + newPolyLine = True + else: + xyz.append([float(val) for val in line.split()]) + nvert = nvert + 1 + if not newPolyLine: + segments.append([nvert - 1, nvert]) + newPolyLine = False + +xyz = np.asarray(xyz) +# add extra column for z coordinates +xyz = np.insert(xyz, xyz.shape[1], 0, axis=1) + +segments = np.asarray(segments) + +if args.proj: + from pyproj import Transformer + + f = lambda x: {"proj": "geocent", "ellps": "WGS84", "datum": "WGS84"} if x == "geocent" else x + transformer = Transformer.from_crs("epsg:4326", f(args.proj[0]), always_xy=True) + xyz[:, 0], xyz[:, 1], xyz[:, 2] = transformer.transform(xyz[:, 0], xyz[:, 1], xyz[:, 2]) + if args.proj[0] != "geocent": + xyz[:, 2] = args.z[0] +else: + xyz[:, 2] = args.z[0] + +xyz[:, 0] -= args.recenter[0] +xyz[:, 1] -= args.recenter[1] + +nvert = xyz.shape[0] +nseg = segments.shape[0] + + +# Now write vtk file +with open("CoastLine.vtk", "w") as fout: + fout.write( + f"""# vtk DataFile Version 2.0 +parabola - polyline +ASCII +DATASET POLYDATA +POINTS {nvert} float +""" + ) + np.savetxt(fout, xyz, "%e %e %e") + fout.write(f"\nLINES {nseg} {3*nseg}\n") + np.savetxt(fout, segments - 1, "2 %d %d") + fout.write("\n") + +print("CoastLine.vtk successfully created") diff --git a/tools/Create_movie.m b/tools/Create_movie.m new file mode 100644 index 0000000..ccf2a3c --- /dev/null +++ b/tools/Create_movie.m @@ -0,0 +1,86 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2009, SeisSol Group +% All rights reserved. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are met: +% +% 1. Redistributions of source code must retain the above copyright notice, +% this list of conditions and the following disclaimer. +% +% 2. Redistributions in binary form must reproduce the above copyright notice, +% this list of conditions and the following disclaimer in the documentation +% and/or other materials provided with the distribution. +% +% 3. Neither the name of the copyright holder nor the names of its +% contributors may be used to endorse or promote products derived from this +% software without specific prior written permission. +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. +% +% @section DESCRIPTION +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% Create_movie builds avi-movies from a list of numbered %%') +disp(' %% images (e.g. jpg). %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' '),disp(' ') +clear, close all; + +%input data +filename = input(' Specify root-filename : ','s'); +fps = input(' Specify frames per sec : '); +nf = input(' Specify number of frames: '); + +fig = figure; +M = moviein(nf); +set(fig,'DoubleBuffer','on'); +mov = avifile([filename,'.avi'],'fps',fps,'compression','none','quality',100); + +counter = 0; + +%start loop over images +for i = 1:nf + + counter = counter + 1; + + eval(['tmp = imread(''',filename,num2str(i),''',''jpg'');']); + + h = image(tmp); + + %Annotation of text if necessary (e.g. frames show every 5 sec starting form 0 sec) + %text(300,500,['Time: ',num2str(5*(i-1)),'sec']) + + drawnow; + + M(:,counter)=getframe; + F = getframe(gca); + mov = addframe(mov,F); + +end + +mov = close(mov); + +disp('-----------------------------------------------------------------------------------') +disp(sprintf('\t Movie built successfully !')); +disp('-----------------------------------------------------------------------------------') +disp(' ') diff --git a/tools/Plot_seissol_snapshot.m b/tools/Plot_seissol_snapshot.m new file mode 100644 index 0000000..a0904ce --- /dev/null +++ b/tools/Plot_seissol_snapshot.m @@ -0,0 +1,92 @@ +%% +% @file +% This file is part of SeisSol. +% +% @author Martin Kaeser (martin.kaeser AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/kaeser) +% +% @section LICENSE +% Copyright (c) 2008, SeisSol Group +% All rights reserved. +% +% Redistribution and use in source and binary forms, with or without +% modification, are permitted provided that the following conditions are met: +% +% 1. Redistributions of source code must retain the above copyright notice, +% this list of conditions and the following disclaimer. +% +% 2. Redistributions in binary form must reproduce the above copyright notice, +% this list of conditions and the following disclaimer in the documentation +% and/or other materials provided with the distribution. +% +% 3. Neither the name of the copyright holder nor the names of its +% contributors may be used to endorse or promote products derived from this +% software without specific prior written permission. +% +% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +% POSSIBILITY OF SUCH DAMAGE. +% +% @section DESCRIPTION + +home; +disp(' ') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' %% %%') +disp(' %% Plot_SEISSOL_Snapshots reads in the TECPLOT format %%') +disp(' %% of the typical SEISSOL snapshot output and visualizes %%') +disp(' %% it on a color-coded triangulation surface plot. %%') +disp(' %% %%') +disp(' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%') +disp(' '),disp(' ') + +clear, close all; +filename = input(' Specify root-filename: ','s'); +nsd = input(' Give number of subdomain: '); +var = input(' Specify variable number (1-5): '); +disp(' '), +disp('ATTENTION: Make sure to have 10 data columns in your .tri. file or change it in this matlab script!') +disp(' ') +disp(' '), disp(' Creating seismogram-processor relation ...' ) + +for i=1:nsd + + %define character string of input files + if(i<=10) + in_file = [filename,'.000',num2str(i-1),'.tri.dat']; + elseif(i>10 & i<=100) + in_file = [filename,'.00',num2str(i-1),'.tri.dat']; + else + in_file = [filename,'.0',num2str(i-1),'.tri.dat']; + end + + %open file and read data + fid = fopen(in_file); + + junk = fgetl(fid); time = str2num(junk(24:end-1)); + junk = fgetl(fid); + junk = fgetl(fid); nX = str2num(junk(9:20)); nE = str2num(junk(26:37)); + data = fscanf(fid,'%g',[10,nX]); data = data'; % change here if the .tri. file does not contain 10 columns! + %extract vertices + X = data(:,1:2); data(:,1:2) = []; + %extrace connectivity matrix + con_mat = fscanf(fid,'%g',[3,nE]); con_mat = con_mat'; + %visualize results + trisurf(con_mat,X(:,1),X(:,2),data(:,var)); hold on; + + disp(sprintf(' Visualizing subdomain %d',i)); + + fclose(fid); + + xlabel('x [m]'), ylabel('y [m]'), zlabel(['Variable ',num2str(var)]); view(0,90) + +end + +disp(sprintf(' \n Finished! \n')); diff --git a/tools/animate_2d_output.py b/tools/animate_2d_output.py new file mode 100644 index 0000000..c55446d --- /dev/null +++ b/tools/animate_2d_output.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import animation +from seissolxdmf import seissolxdmf as sx + +plt.rcParams['savefig.facecolor'] = 'w' + + +def main(): + description = '''Creates an animation of the velocity/displacement + wavefield from SeisSol free-surface XDMF output file.''' + parser = argparse.ArgumentParser(description=description) + parser.add_argument( + 'fin', help='SeisSol surface (XDMF) output file.') + parser.add_argument( + 'fout', help='Output animation file.', default='animation.gif') + parser.add_argument( + 'dataid', help='Which data component to use', default='v1') + parser.add_argument( + '--cmap', help='Matplotlib colormap', default='bwr') + parser.add_argument( + '--dpi', default=200, type=int) + parser.add_argument( + '--bounds', nargs='+', help='xmin xmax ymin ymax', type=float) + args = parser.parse_args() + + ds = sx(args.fin) + data = ds.ReadData(args.dataid) + geo = ds.ReadGeometry() + connect = ds.ReadConnect() + dt = ds.ReadTimeStep() + nt = data.shape[0] + tarr = np.arange(0.0, nt * dt, dt) + + fig = plt.figure(figsize=(8, 7)) + ax = plt.subplot(111) + dmax = abs(data).max() + im = ax.tripcolor( + geo.T[0], geo.T[1], connect, data[0], cmap=args.cmap, vmin=-dmax, + vmax=dmax) + ax.set_xlabel('x (m)') + ax.set_ylabel('y (m)') + ax.set_xlim(args.bounds[0], args.bounds[1]) + ax.set_ylim(args.bounds[2], args.bounds[3]) + fig.colorbar(im, label=args.dataid) + + def animate(i): + im.set_array(data[i]) + ax.set_title('$t=$%.1f' % tarr[i]) + + def progress_callback(i, n): + return print(f'Saving frame {i} of {n}') + + anim = animation.FuncAnimation( + fig, animate, frames=data.shape[0], + interval=dt * 1000) + anim.save(args.fout, dpi=args.dpi, progress_callback=progress_callback) + plt.close('all') + + +if __name__ == '__main__': + main() diff --git a/tools/extractDataFromUnstructuredOutput.py b/tools/extractDataFromUnstructuredOutput.py new file mode 100644 index 0000000..5fd172c --- /dev/null +++ b/tools/extractDataFromUnstructuredOutput.py @@ -0,0 +1,219 @@ +import os +import os.path +import seissolxdmf +import shutil +import recreateXdmf +import numpy as np +import argparse + +parser = argparse.ArgumentParser(description="resample output file and write as binary files") +parser.add_argument("xdmfFilename", help="xdmf output file") +parser.add_argument("--add2prefix", help="string to append to prefix for new file", type=str, default="_resampled") +parser.add_argument("--Data", nargs="+", metavar=("variable"), help="Data to resample (example SRs, or all)", required=True) +parser.add_argument("--downsample", help="write one out of n output", type=int) +parser.add_argument("--precision", type=str, choices=["float", "double"], default="float", help="precision of output file") +parser.add_argument("--backend", type=str, choices=["hdf5", "raw"], default="hdf5", help="backend used: raw (.bin file), hdf5 (.h5)") +parser.add_argument("--last", dest="last", default=False, action="store_true", help="output last time step") +parser.add_argument("--idt", nargs="+", help="list of time step to write (ex $(seq 7 3 28))", type=int) +parser.add_argument( + "--xfilter", + nargs=2, + metavar=("xmin", "xmax"), + help="output only cells with x center coordinate in range xmin xmax", + type=float, +) +parser.add_argument( + "--yfilter", + nargs=2, + metavar=("ymin", "ymax"), + help="output only cells with y center coordinate in range ymin ymax", + type=float, +) +parser.add_argument( + "--zfilter", + nargs=2, + metavar=("zmin", "zmax"), + help="output only cells with z center coordinate in range zmin zmax", + type=float, +) +args = parser.parse_args() + + +class seissolxdmfExtended(seissolxdmf.seissolxdmf): + def ReadData(self, dataName, idt=-1): + if dataName == "SR": + SRs = super().ReadData("SRs", idt) + SRd = super().ReadData("SRd", idt) + return np.sqrt(SRs ** 2 + SRd ** 2) + else: + return super().ReadData(dataName, idt) + + +sx = seissolxdmfExtended(args.xdmfFilename) +xyz = sx.ReadGeometry() +connect = sx.ReadConnect() + +spatial_filtering = (args.xfilter or args.yfilter) or args.zfilter +nElements = connect.shape[0] + +if spatial_filtering: + print("Warning: spatial filtering significantly slows down this script") + ids = range(0, sx.nElements) + xyzc = (xyz[connect[:, 0], :] + xyz[connect[:, 1], :] + xyz[connect[:, 2], :]) / 3.0 + + def filter_cells(coords, filter_range): + m = 0.5 * (filter_range[0] + filter_range[1]) + d = 0.5 * (filter_range[1] - filter_range[0]) + return np.where(np.abs(coords[:] - m) < d)[0] + + if args.xfilter: + id0 = filter_cells(xyzc[:, 0], args.xfilter) + ids = np.intersect1d(ids, id0) if len(ids) else id0 + if args.yfilter: + id0 = filter_cells(xyzc[:, 1], args.yfilter) + ids = np.intersect1d(ids, id0) if len(ids) else id0 + if args.zfilter: + id0 = filter_cells(xyzc[:, 2], args.zfilter) + ids = np.intersect1d(ids, id0) if len(ids) else id0 + + if len(ids): + connect = connect[ids, :] + nElements = connect.shape[0] + if nElements != sx.nElements: + print(f"extracting {nElements} cells out of {sx.nElements}") + else: + spatial_filtering = False + else: + raise ValueError("all elements are outside filter range") +ndt = sx.ndt + +if args.last: + indices = [ndt - 1] + if args.idt or args.downsample: + print("last option cannot be used together with idt and downsample options") + exit() +else: + if args.idt and args.downsample: + print("idt and downsample options cannot be used together") + exit() + elif not args.downsample: + indices = args.idt + else: + indices = range(0, ndt, args.downsample) + +# Check if input is in hdf5 format or not +first_data_field = list(sx.ReadAvailableDataFields())[0] +dataLocation, data_prec, MemDimension = sx.GetDataLocationPrecisionMemDimension(first_data_field) +splitArgs = dataLocation.split(":") +if len(splitArgs) == 2: + isHdf5 = True +else: + isHdf5 = False + +if args.precision == "double": + myDtype = "float64" + myprec = 8 +else: + myDtype = "float32" + myprec = 4 + +if args.backend == "raw": + write2Binary = True +else: + write2Binary = False + +prefix = os.path.splitext(args.xdmfFilename)[0] +prefix_new = recreateXdmf.generate_new_prefix(prefix, args.add2prefix) + +# Create folders if necessary +if write2Binary: + if not os.path.exists(prefix_new + "_cell"): + os.mkdir(prefix_new + "_cell") + if not os.path.exists(prefix_new + "_cell/mesh0/"): + os.mkdir(prefix_new + "_cell/mesh0/") + if not os.path.exists(prefix_new + "_vertex"): + os.mkdir(prefix_new + "_vertex") + if not os.path.exists(prefix_new + "_vertex/mesh0/"): + os.mkdir(prefix_new + "_vertex/mesh0/") + + +# Write geometry and connect +if write2Binary: + fn2 = prefix_new + "_cell/mesh0/connect.bin" + if isHdf5: + # write Connect + output_file = open(fn2, "wb") + connect.tofile(output_file) + output_file.close() + else: + shutil.copy2(os.path.splitext(args.xdmfFilename)[0] + "_cell/mesh0/connect.bin", fn2) + print("done writing " + fn2) + # write geometry + fn3 = prefix_new + "_vertex/mesh0/geometry.bin" + output_file = open(fn3, "wb") + xyz.tofile(output_file) + output_file.close() + print("done writing " + fn3) +else: + import h5py + + # write geometry to hdf5 format + h5fv = h5py.File(prefix_new + "_vertex.h5", "w") + h5fv.create_dataset("/mesh0/geometry", data=xyz) + h5fv.close() + print("done writing " + prefix_new + "_vertex.h5") + # write connect to hdf5 format + h5fc = h5py.File(prefix_new + "_cell.h5", "w") + h5fc.create_dataset("/mesh0/connect", data=connect) + + +# Write data items +if args.Data[0]=='all': + args.Data = sorted(sx.ReadAvailableDataFields()) + for to_remove in ["partition", "locationFlag"]: + if to_remove in args.Data: + args.Data.remove(to_remove) + print(f"args.Data was set to all and now contains {args.Data}") + +for ida, sdata in enumerate(args.Data): + if write2Binary: + fname2 = prefix_new + "_cell/mesh0/" + args.Data[ida] + ".bin" + output_file = open(fname2, "wb") + else: + dset = h5fc.create_dataset("/mesh0/" + args.Data[ida], (len(indices), nElements), dtype=myDtype) + # read only one row + print(sdata, end=" ", flush=True) + for kk, i in enumerate(indices): + if (kk % 10 == 0) and kk > 0: + print(kk) + else: + print(kk, end=" ", flush=True) + if i >= ndt: + print("ignoring index %d>=ndt=%d" % (i, ndt)) + continue + if spatial_filtering: + myData = sx.ReadData(args.Data[ida], idt=i)[ids] + else: + myData = sx.ReadData(args.Data[ida], idt=i) + if write2Binary: + myData.astype(myDtype).tofile(output_file) + else: + dset[kk, :] = myData[:] + if write2Binary: + output_file.close() + print("done writing " + fname2) + +if not write2Binary: + h5fc.close() + print("done writing " + prefix_new + "_cell.h5") + +# Now recreate the Xdmf + +prefix = os.path.splitext(args.xdmfFilename)[0] + +# Read all parameters from the xdmf file of SeisSol (if any) + +dt = recreateXdmf.ReadDtFromXdmf(args.xdmfFilename) +nvertex = recreateXdmf.ReadNvertexFromXdmf(args.xdmfFilename) +ndt, nmem = recreateXdmf.ReadNdtNmemFromXdmf(args.xdmfFilename) +recreateXdmf.recreateXdmf(prefix, prefix_new, nvertex, nElements, nElements, dt, indices, args.Data, not write2Binary, myprec, args.add2prefix) diff --git a/tools/generate_rec_vtp.py b/tools/generate_rec_vtp.py new file mode 100644 index 0000000..533a26c --- /dev/null +++ b/tools/generate_rec_vtp.py @@ -0,0 +1,108 @@ +#! /usr/bin/python +## +# @file +# This file is part of SeisSol. +# +# @author Stefan Wenk (wenk AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/wenk) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# Creates vtk file for receiver coordinates from SeisSol parameter file +# +# usage: +# python generate_rec_vtk +# i.e. python generate_rec_vtk.py PARAMETER.par +# +# output: +# rec.vtp receiver coordinates in vtk format + +import sys + +try: + fi = open(sys.argv[1], 'r') +except: + sys.exit('please give SeisSol parameter file and call: python generate_rec_vtk.py PARAMETER.par') + +count = 0 +for line in fi: + if line[0] == '<': + count = count+1 + if count == 27: + break + +for i,line in enumerate(fi): + if i == 10: + nrec = int(line.split()[0]) + break + +a = [] +for i,line in enumerate(fi): + a = a+[line] + if i == nrec-1: + break + +fo = open('rec.vtp', 'w') + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +for i in a: + fo.write(i) + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +for i in range(nrec): + fo.write(str(i)+' ') + +fo.write('\n') +fo.write('') +fo.write('\n') +for i in range(1,nrec+1): + fo.write(str(i)+' ') + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +for i in range(1,nrec+1): + fo.write(str(i)+'\n') + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') diff --git a/tools/neu2vtk.py b/tools/neu2vtk.py new file mode 100644 index 0000000..0eabf1e --- /dev/null +++ b/tools/neu2vtk.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Martin van Driel (driel AT geophysik.uni-muenchen.de, http://www.geophysik.lmu.de/Members/driel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# neu2vtk +# +# creates vtk-file 'neufile.vtk' +# +# for tetrahedral meshes only!!!! + +import sys + +fi = open(sys.argv[1], 'r') +fo = open(sys.argv[1] + '.vtk', 'w') + +for i in range(7): + s = fi.readline() +s = s.split() + +npts = int(s[0]) +nelem = int(s[1]) + +fo.write('# vtk DataFile Version 3.0\n') +fo.write('vtk ' + sys.argv[1] + '\n') +fo.write('ASCII\n') +fo.write('DATASET UNSTRUCTURED_GRID\n') + +fo.write('POINTS ' + str(npts) + ' float\n') + +fi.readline() +fi.readline() + +for i in range(npts): + s = fi.readline().split() + fo.write(s[1] + ' ' + s[2] + ' ' + s[3] + '\n') + +fi.readline() +fi.readline() + +fo.write('CELLS ' + str(nelem) + ' ' + str(nelem*5) + '\n') + +for i in range(nelem): + s = fi.readline().split() + s3 = str(int(s[3]) - 1) + s4 = str(int(s[4]) - 1) + s5 = str(int(s[5]) - 1) + s6 = str(int(s[6]) - 1) + fo.write('4 ' + s3 + ' ' + s4 + ' ' + s5 + ' ' + s6 + '\n') + +fo.write('CELL_TYPES ' + str(nelem) + '\n') + +for i in range(nelem): + fo.write('10\n') + +fi.close() +fo.close() diff --git a/tools/pickpoint2mseed.py b/tools/pickpoint2mseed.py new file mode 100755 index 0000000..d6c6911 --- /dev/null +++ b/tools/pickpoint2mseed.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Martin van Driel (driel AT geophysik.uni-muenchen.de, http://www.geophysik.lmu.de/Members/driel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python pickpoint2mseed +# +# creates two files: +# .coords contains the station coordinates +# .mseed contains the seismograms in mini-seed format +# +# requires obspy (http://www.obspy.org) + +import numpy as np +import sys +import glob +try: + from obspy.core import Trace, Stream, UTCDateTime +except: + sys.exit('please install obspy (http://www.obspy.org)') + +t = UTCDateTime(0) + +try: + files = sys.argv[1] +except: + sys.exit('usage: python pickpoint2mseed ') + +traces = [] + +f = open(files + '.coords', 'w') + +for file in glob.iglob(files + '*.dat'): + print file + header = [x.replace("\"","").replace(","," ").replace("="," ").split() for x in open(file).readlines()[:5]] + stationID = header[0][-1] + channels = header[1][2:] + x1 = header[2][2] + x2 = header[3][2] + x3 = header[4][2] + + if x1[0] != '-': + x1 = ' ' + x1 + if x2[0] != '-': + x2 = ' ' + x2 + if x3[0] != '-': + x3 = ' ' + x3 + + f.write(stationID[-4:] + ' ' + x1 + ' ' + x2 + ' ' + x3 + '\n') + + dat = np.loadtxt(file, skiprows=5) + + for i, chan in enumerate(channels): + stats = {'network': 'SG', + 'station': stationID[-4:], + 'location': '', + 'channel': chan, + 'npts': len(dat[:,i+1]), + 'sampling_rate': 1./(dat[1,0] - dat[0,0]), + 'starttime': t, + 'mseed' : {'dataquality': 'D'}} + traces.append(Trace(data=dat[:,i+1], header=stats)) + + +f.close() + +st = Stream(traces) +st.sort() +print st +fname = files + '.mseed' +print fname +st.write(fname, format='MSEED') diff --git a/tools/plot_imt_map.py b/tools/plot_imt_map.py new file mode 100644 index 0000000..5d2cf57 --- /dev/null +++ b/tools/plot_imt_map.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python3 +import argparse +import numpy as np +import matplotlib.pyplot as plt +from seissolxdmf import seissolxdmf as sx + + +def main(): + description = '''Plots a map view of the ground motion intensity.''' + parser = argparse.ArgumentParser(description=description) + parser.add_argument( + 'fin', help='XDMF input file name.') + parser.add_argument( + 'fout', help='Output image file name.') + parser.add_argument( + 'dataid', help='Which intensity metric to use', default='PGV') + parser.add_argument( + '--cmap', help='Matplotlib colormap', default='viridis') + parser.add_argument( + '--dpi', default=200) + parser.add_argument( + '--log', help='Plot log of intensity', action='store_true', + default=False) + parser.add_argument( + '--bounds', nargs='+', help='xmin xmax ymin ymax', type=float) + parser.add_argument( + '--contour', action='store_true', default=False) + args = parser.parse_args() + + ds = sx(args.fin) + geo = ds.ReadGeometry() + connect = ds.ReadConnect() + coords = geo[connect].mean(axis=1) + x, y = coords.T[0], coords.T[1] + imt = ds.ReadData(args.dataid) + clabel = args.dataid + if args.log: + imt = np.log(imt) + clabel = 'log(%s)' % clabel + + if args.bounds: + idx = ((x > args.bounds[0]) & (x < args.bounds[1]) & + (y > args.bounds[2]) & (y < args.bounds[3])) + vmin = imt[idx].min() + vmax = imt[idx].max() + else: + vmin, vmax = None, None + + plt.figure(figsize=(6, 5)) + + if args.contour: + plt.tricontourf(x, y, imt, vmin=vmin, vmax=vmax) + else: + plt.tripcolor(geo.T[0], geo.T[1], connect, imt, vmin=vmin, vmax=vmax) + if args.bounds: + plt.xlim(args.bounds[0], args.bounds[1]) + plt.ylim(args.bounds[2], args.bounds[3]) + plt.colorbar(label=clabel) + plt.xlabel('x (m)') + plt.ylabel('y (m)') + plt.savefig(args.fout, dpi=int(args.dpi)) + print('Image file saved to %s' % args.fout) + plt.close('all') + + +if __name__ == '__main__': + main() diff --git a/tools/recreateXdmf.py b/tools/recreateXdmf.py new file mode 100644 index 0000000..de3e816 --- /dev/null +++ b/tools/recreateXdmf.py @@ -0,0 +1,193 @@ +import argparse +import subprocess +import xml.etree.ElementTree +import os + + +def next_power_of_2(x): + return 1 if x == 0 else 2 ** ((x - 1).bit_length()) + + +def generate_new_prefix(prefix, append2prefix): + if append2prefix == "": + append2prefix = "_resampled" + prefix = os.path.basename(prefix) + lsplit = prefix.split("-") + if len(lsplit) > 1: + if lsplit[-1] in ["surface", "low", "fault"]: + prefix0 = "-".join(lsplit[0:-1]) + prefix_new = prefix0 + append2prefix + "-" + lsplit[-1] + else: + prefix_new = prefix + append2prefix + return prefix_new + + +def recreateXdmf(prefix, prefix_new, nvertex, ncells, nmem, dt, indices, lsData, tohdf5=False, prec=8, append2prefix="_resampled"): + full_prefix = prefix + prefix = os.path.basename(prefix) + prefix_new = os.path.basename(prefix_new) + + # fault and surface output have 3 colums in connect + ncolConnect = 4 + scell = "Tetrahedron" + lsplit = prefix.split("-") + if len(lsplit) > 1: + if lsplit[-1] in ["surface", "fault"]: + ncolConnect = 3 + scell = "Triangle" + + if tohdf5: + colonOrNothing = ".h5:" + DataExtension = "" + DataFormat = "HDF" + else: + colonOrNothing = "" + DataExtension = ".bin" + DataFormat = "Binary" + + # create and print the new Xdmf file + header = """ + + + + """ + + score = "" + + for i, ii in enumerate(indices): + score = ( + score + + """ + + %s_cell%s/mesh0/connect%s + + + %s_vertex%s/mesh0/geometry%s + + \n" + + score = ( + score + + """ + +""" + ) + prefix0 = generate_new_prefix(full_prefix, append2prefix) + fout = open(prefix0 + ".xdmf", "w") + fout.write(header) + fout.write("\n") + fout.write(score) + fout.write("\n") + fout.close() + print("done writing " + prefix0 + ".xdmf") + + +def ReadNcellsFromXdmf(xdmfFile): + out = subprocess.check_output(["grep connect " + xdmfFile + " | head -n1"], shell=True) + e = xml.etree.ElementTree.fromstring(out) + dimstring = e.attrib["Dimensions"].split() + return int(dimstring[0]) + + +def ReadDtFromXdmf(xdmfFile): + out = subprocess.check_output(["grep Time " + xdmfFile + " | head -n3| tail -n1"], shell=True) + e = xml.etree.ElementTree.fromstring(out) + dimstring = e.attrib["Value"].split() + return float(dimstring[0]) + + +def ReadNvertexFromXdmf(xdmfFile): + out = subprocess.check_output(["grep geometry " + xdmfFile + " | head -n1"], shell=True) + e = xml.etree.ElementTree.fromstring(out) + dimstring = e.attrib["Dimensions"].split() + return int(dimstring[0]) + + +def ReadNdtNmemFromXdmf(xdmfFile): + out = subprocess.check_output(["grep DataItem " + xdmfFile + " | tail -n2 | head -n1"], shell=True) + e = xml.etree.ElementTree.fromstring(out) + dimstring = e.attrib["Dimensions"].split() + # return (int(dimstring[0])-1, int(dimstring[1])) + return (int(dimstring[0]), int(dimstring[1])) + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser(description="recreate a xdmf file") + parser.add_argument("prefix", help="prefix including -fault or -surface") + parser.add_argument("--idt", nargs="+", help="list of time step to differenciate (1st = 0); -1 = all", type=int) + parser.add_argument("--nvertex", nargs=1, metavar=("nvertex"), help="number of vertex (read if not given)", type=int) + parser.add_argument("--ncells", nargs=1, metavar=("ncells"), help="number of cells (read if not given)", type=int) + parser.add_argument("--dt", nargs=1, metavar=("dt"), help="output time step (read if not given)", type=float) + parser.add_argument("--ndt", nargs=1, metavar=("ndt"), help="number of time steps to output (read if not given)", type=int) + parser.add_argument("--Data", nargs="+", help="list of data variable to write") + args = parser.parse_args() + + prefix = args.prefix + xdmfFile = prefix + ".xdmf" + ### Read all parameters from the xdmf file of SeisSol (if any) + if args.ncells != None: + ncells = args.ncells[0] + nmem = next_power_of_2(ncells) + # nmem = next_power_of_2(ncells)*2 + else: + ncells = ReadNcellsFromXdmf(xdmfFile) + + if args.dt != None: + dt = args.dt[0] + else: + dt = ReadDtFromXdmf(xdmfFile) + + if args.nvertex != None: + nvertex = args.nvertex[0] + else: + nvertex = ReadNvertexFromXdmf(xdmfFile) + + if args.ndt != None: + ndt = args.ndt[0] + else: + ndt, nmem = ReadNdtNmemFromXdmf(xdmfFile) + + recreateXdmf(prefix, prefix, nvertex, ncells, nmem, dt, range(ndt), args.Data) diff --git a/tools/sismowine2receiver.py b/tools/sismowine2receiver.py new file mode 100644 index 0000000..6b0d7db --- /dev/null +++ b/tools/sismowine2receiver.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Carsten Uphoff (c.uphoff AT tum.de, http://www5.in.tum.de/wiki/index.php/Carsten_Uphoff,_M.Sc.) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# sismowine2receiver + +import sys +import numpy + +class SismowineFile: + def __init__(self, filename): + f = open(filename, 'r') + header = f.readline() + header = header.split() + self.nsamples = int(header[0]) + self.dt = float(header[1]) + self.nrec = int(header[2]) + self.data = numpy.empty([self.nsamples, self.nrec], order='F', dtype=float) + for rec in range(self.nrec): + for sample in range(self.nsamples): + self.data[sample, rec] = float(f.readline()) + + + +xdat = SismowineFile(sys.argv[1] + '/x.dat') +ydat = SismowineFile(sys.argv[1] + '/y.dat') +zdat = SismowineFile(sys.argv[1] + '/z.dat') + +for rec in range(xdat.nrec): + fo = open('{0}-{1:03d}.dat'.format(sys.argv[2], rec+1), 'w') + fo.write('TITLE = "Reference for receiver number {0:03d}"\n'.format(rec)) + fo.write('VARIABLES = "Time","u","v","w"\n') + fo.write('# x1 NaN\n') + fo.write('# x2 NaN\n') + fo.write('# x3 NaN\n') + for sample in range(xdat.nsamples): + fo.write('{} {} {} {}\n'.format(sample * xdat.dt, xdat.data[sample, rec], ydat.data[sample, rec], zdat.data[sample, rec])) diff --git a/tools/tecp2vtk.py b/tools/tecp2vtk.py new file mode 100644 index 0000000..067ebea --- /dev/null +++ b/tools/tecp2vtk.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Martin van Driel (driel AT geophysik.uni-muenchen.de, http://www.geophysik.lmu.de/Members/driel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python tecp2vtk +# +# creates vtk-file 'tecpfile.vtk' +# +# for tetrahedral meshes only!!!! + +import sys +import numpy as np + +fi = open(sys.argv[1], 'r') +fo = open(sys.argv[1] + '.vtk', 'w') + +title = fi.readline().split('"')[1].strip() + +vars = fi.readline().split()[2:] + +for i, var in enumerate(vars): + vars[i] = var.strip('"') + +nvar = len(vars) +s = fi.readline().split() +npts = int(s[2]) +nelem = int(s[4]) + +data = np.empty((nvar, npts)) + +for i in range(npts): + data[:,i] = np.array(fi.readline().split(), dtype=float) + +elems = np.zeros((4, nelem), dtype=int) +for i in range(nelem): + elems[:,i] = np.array(fi.readline().split(), dtype=int) - 1 + +fi.close() + + +fo.write('# vtk DataFile Version 3.0\n') +fo.write(title + '\n') +fo.write('ASCII\n') +fo.write('DATASET UNSTRUCTURED_GRID\n') + +fo.write('POINTS ' + str(npts) + ' float\n') + + +for i in range(npts): + fo.write(str(data[0,i]) + ' ' + str(data[1,i]) + ' ' + str(data[2,i]) + '\n') + +fo.write('CELLS ' + str(nelem) + ' ' + str(nelem*5) + '\n') + +for i in range(nelem): + fo.write('4 ' + str(elems[0,i]) + ' ' + str(elems[1,i]) + ' ' + str(elems[2,i]) + ' ' + str(elems[3,i]) + '\n') + +fo.write('CELL_TYPES ' + str(nelem) + '\n') + +for i in range(nelem): + fo.write('10\n') + +fo.write('POINT_DATA ' + str(npts) + '\n') +for j, var in enumerate(vars[3:]): + fo.write('SCALARS ' + var + ' float 1\n') + fo.write('LOOKUP_TABLE default\n') + + for i in range(npts): + fo.write(str(data[j+3,i]) + '\n') + +fo.close() diff --git a/tools/tecp2vtk_merge.py b/tools/tecp2vtk_merge.py new file mode 100644 index 0000000..f9b0c78 --- /dev/null +++ b/tools/tecp2vtk_merge.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Martin van Driel (driel AT geophysik.uni-muenchen.de, http://www.geophysik.lmu.de/Members/driel) +# @author Alice Gabriel (gabriel AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/gabriel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python tecp2vtk_merge [-b] +# i.e. python tecp2vtk_merge.py 000000126 +# +# creates vtk-file '.vtk' +# +# options: +# -b: binary output +# +# for tetrahedral meshes only!!! +# +# requires pyvtk (http://pypi.python.org/pypi/PyVTK) + +import sys +import numpy as np +import glob +import os + +try: + import pyvtk +except: + sys.exit('please install pyvtk (http://pypi.python.org/pypi/PyVTK)') + +binary = False +tetfiles = [] + +if len(sys.argv) == 2: + files = sys.argv[1] +elif len(sys.argv) == 3: + if sys.argv[1] == '-b': + binary = True + files = sys.argv[2] +else: + sys.exit('usage: python tecp2vtk_merge [-b] ') + +for path, dirs, fls in os.walk('./'): + for f in fls: + if f.endswith(files + '.tet.dat'): + tetfiles.append(os.path.join(path, f)) + +if len(tetfiles) == 0: + sys.exit('no such file: timestep-' + files + '.tet.dat') + +elems = [] +datal = [] +ntelem = 0 +ntpts = 0 + +for file in tetfiles: + print 'reading ' + file + fi = open(file, 'r') + + title = fi.readline().split('"')[1].strip() + + vars = fi.readline().split()[2:] + + for i, var in enumerate(vars): + vars[i] = var.strip('"') + + nvar = len(vars) + s = fi.readline().split() + npts = int(s[2]) + nelem = int(s[4]) + + for i in range(npts): + datal.append(fi.readline().split()) + + for i in range(nelem): + elems.append((np.array(fi.readline().split(), dtype=int) - 1 + ntpts).tolist()) + + ntelem += nelem + ntpts += npts + + fi.close() + +data = np.array(datal, dtype=float).T +elems = np.array(elems, dtype=int).T + +pl = [] +for i in np.arange(data.shape[1]): + pl.append(tuple(data[:3,i])) + +el = [] +for i in np.arange(elems.shape[1]): + el.append(tuple(elems[:,i])) + +grid = pyvtk.UnstructuredGrid(pl, tetra=el) + +pointdata = pyvtk.PointData() +for j, var in enumerate(vars[3:]): + pointdata.append(pyvtk.Scalars(data[j+3,:].tolist(), name=var, lookup_table='default')) + +vtk = pyvtk.VtkData(grid, pointdata, title) +if binary: + vtk.tofile(files + '.vtk', 'binary') +else: + vtk.tofile(files + '.vtk', 'ascii') + +print files + '.vtk written' diff --git a/tools/tecp2vtk_merge_loop.py b/tools/tecp2vtk_merge_loop.py new file mode 100644 index 0000000..8c22129 --- /dev/null +++ b/tools/tecp2vtk_merge_loop.py @@ -0,0 +1,155 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Martin van Driel (driel AT geophysik.uni-muenchen.de, http://www.geophysik.lmu.de/Members/driel) +# @author Simon Kremers (kremers AT geophysik.uni-muenchen.de) +# @author Alice Gabriel (gabriel AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/gabriel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python tecp2vtk_merge [-b] +# i.e. python tecp2vtk_merge_loop.py data +# +# creates vtk-files '.vtk' for all timesteps in subfolders +# +# options: +# -b: binary output +# +# for tetrahedral meshes only!!! +# +# requires pyvtk (http://pypi.python.org/pypi/PyVTK) + +import sys +import numpy as np +import glob +import os + +try: + import pyvtk +except: + sys.exit('please install pyvtk (binary = Falsehttp://pypi.python.org/pypi/PyVTK)') + +binary = False +tetfiles = [] +ts = [] + +if len(sys.argv) == 2: + files = sys.argv[1] +elif len(sys.argv) == 3: + if sys.argv[1] == '-b': + binary = True + files = sys.argv[2] +else: + sys.exit('usage: python tecp2vtk_merge [-b] (e.g. "data" )') + +for path, dirs, fls in os.walk('./'): + for d in dirs: + if d.startswith(files): + for f in glob.iglob(os.path.join(path, d, '*.tet.dat')): + tetfiles.append(f) + + +if len(tetfiles) == 0: + sys.exit('no such file: ./' + files + '*/*.tet.dat') + +# get number and name of timesteps +for path, dirs, fls in os.walk('./'): + for d in dirs: + if d.startswith(files): + for f in glob.iglob(os.path.join(path, d, '*.tet.dat')): + ts.append(f) + print ts + break + +# loop over timesteps +for time in ts: + t_step = (time.split('-')[2]).split('.')[0] + elems = [] + datal = [] + ntelem = 0 + ntpts = 0 + + for file in tetfiles: + if file.endswith(t_step + '.tet.dat'): + print 'reading ' + file + fi = open(file, 'r') + + title = fi.readline().split('"')[1].strip() + + vars = fi.readline().split()[2:] + + for i, var in enumerate(vars): + vars[i] = var.strip('"') + + nvar = len(vars) + s = fi.readline().split() + npts = int(s[2]) + nelem = int(s[4]) + + for i in range(npts): + datal.append(fi.readline().split()) + + for i in range(nelem): + elems.append((np.array(fi.readline().split(), dtype=int) - 1 + ntpts).tolist()) + + ntelem += nelem + ntpts += npts + + fi.close() + + data = np.array(datal, dtype=float).T + elems = np.array(elems, dtype=int).T + + pl = [] + for i in np.arange(data.shape[1]): + pl.append(tuple(data[:3,i])) + + el = [] + for i in np.arange(elems.shape[1]): + el.append(tuple(elems[:,i])) + + grid = pyvtk.UnstructuredGrid(pl, tetra=el) + + pointdata = pyvtk.PointData() + for j, var in enumerate(vars[3:]): + pointdata.append(pyvtk.Scalars(data[j+3,:].tolist(), name=var, lookup_table='default')) + + vtk = pyvtk.VtkData(grid, pointdata, title) + if binary: + vtk.tofile(t_step + '.vtk', 'binary') + else: + vtk.tofile(t_step + '.vtk', 'ascii') + + print t_step + '.vtk written' diff --git a/tools/tecp2vtu.py b/tools/tecp2vtu.py new file mode 100644 index 0000000..6ac75cd --- /dev/null +++ b/tools/tecp2vtu.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Alice Gabriel (gabriel AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/gabriel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python tecp2vtu +# +# creates vtu-file 'tecpfile.vtu' for one tecplot output file +# +# for tetrahedral meshes only! + +import sys +import numpy as np + +fi = open(sys.argv[1], 'r') +fo = open(sys.argv[1] + '.vtu', 'w') + +title = fi.readline().split('"')[1].strip() + +vars = fi.readline().split()[2:] + +for i, var in enumerate(vars): + vars[i] = var.strip('"') + +nvar = len(vars) +s = fi.readline().split() +npts = int(s[2]) +nelem = int(s[4]) + +data = np.empty((nvar, npts)) + +for i in range(npts): + data[:,i] = np.array(fi.readline().split(), dtype=float) + +elems = np.zeros((4, nelem), dtype=int) +for i in range(nelem): + elems[:,i] = np.array(fi.readline().split(), dtype=int) - 1 + +fi.close() + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') + +for i in range(npts): + fo.write(' ' + str('{: .16E}'.format(data[0,i])) + ' ' + str('{: .16E}'.format(data[1,i])) + ' ' + str('{: .16E}'.format(data[2,i])) + '\n') + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') + +for i in range(nelem): + fo.write(' ' + str('{:6d}'.format(elems[0,i])) + ' ' + str('{:6d}'.format(elems[1,i])) + ' ' + str('{:6d}'.format(elems[2,i])) + ' ' + str('{:6d}'.format(elems[3,i])) + '\n') +fo.write('\n') +fo.write('\n') + +for i in range(nelem): + fo.write(' ' + str('{:6d}'.format((i+1)*4)) + '\n') + +fo.write('\n') +fo.write('\n') + +for i in range(nelem): + fo.write(' ' + str('{:6d}'.format(10)) + '\n') + +fo.write('\n') +fo.write('\n') + +fo.write('\n') + +for j, var in enumerate(vars[3:]): + + fo.write('\n') + + for i in range(npts): + fo.write(' ' + str('{: .16E}'.format(data[j+3,i])) + '\n') + + fo.write('\n') + + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.close() diff --git a/tools/tecp2vtu_merge.py b/tools/tecp2vtu_merge.py new file mode 100644 index 0000000..9316f36 --- /dev/null +++ b/tools/tecp2vtu_merge.py @@ -0,0 +1,151 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Alice Gabriel (gabriel AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/gabriel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python tecp2vtk_merge [-b] +# i.e. python tecp2vtk_merge.py data-000000126 +# +# creates vtu-file 'tecpfile.vtu' for one timestep +# +# for tetrahedral meshes only! + +import sys +import numpy as np +import glob + + +fo = open(sys.argv[1] + '.vtu', 'w') + +if len(sys.argv) == 2: + files = sys.argv[1] +else: + sys.exit('usage: python tecp2vtk_merge [-b] ') + +if len(glob.glob(files + '*.tet.dat')) == 0: + sys.exit('no such file: ' + files + '*.tet.dat') + +elems = [] +datal = [] +ntelem = 0 +ntpts = 0 + +for file in glob.iglob(files + '*.tet.dat'): + print 'reading ' + file + fi = open(file, 'r') + + title = fi.readline().split('"')[1].strip() + + vars = fi.readline().split()[2:] + + for i, var in enumerate(vars): + vars[i] = var.strip('"') + + nvar = len(vars) + s = fi.readline().split() + npts = int(s[2]) + nelem = int(s[4]) + + for i in range(npts): + datal.append(fi.readline().split()) + + for i in range(nelem): + elems.append((np.array(fi.readline().split(), dtype=int) - 1 + ntpts).tolist()) + + ntelem += nelem + ntpts += npts + + fi.close() + +data = np.array(datal, dtype=float).T +elems = np.array(elems, dtype=int).T + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') + +for i in range(ntpts): + fo.write(' ' + str('{: .16E}'.format(data[0,i])) + ' ' + str('{: .16E}'.format(data[1,i])) + ' ' + str('{: .16E}'.format(data[2,i])) + '\n') + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') + +for i in range(ntelem): + fo.write(' ' + str('{:6d}'.format(elems[0,i])) + ' ' + str('{:6d}'.format(elems[1,i])) + ' ' + str('{:6d}'.format(elems[2,i])) + ' ' + str('{:6d}'.format(elems[3,i])) + '\n') + +fo.write('') + +fo.write('\n') + +for i in range(ntelem): + fo.write(' ' + str('{:6d}'.format((i+1)*4)) + '\n') + +fo.write('\n') + +fo.write('\n') + + +for i in range(ntelem): + fo.write(' ' + str('{:6d}'.format(10)) + '\n') + +fo.write('\n') + +fo.write('\n') + +fo.write('\n') + + +for j, var in enumerate(vars[3:]): + + fo.write('\n') + + for i in range(npts): + fo.write(' ' + str('{: .16E}'.format(data[j+3,i])) + '\n') + + fo.write('\n') + + +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.write('\n') +fo.close() +print files + '.vtu written' diff --git a/tools/tecp2vtu_merge_loop.py b/tools/tecp2vtu_merge_loop.py new file mode 100644 index 0000000..b9fb01f --- /dev/null +++ b/tools/tecp2vtu_merge_loop.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Alice Gabriel (gabriel AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/gabriel) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# usage: +# python tecp2vtk_merge [-b] +# i.e. python tecp2vtk_merge.py data +# +# creates vtu-files 'tecpfile.vtu' for all timesteps +# +# for tetrahedral meshes only! + +import sys +import numpy as np +import glob + +if len(sys.argv) == 2: + files = sys.argv[1] +else: + sys.exit('usage: python tecp2vtk_merge [-b] ') + +if len(glob.glob(files + '*.tet.dat')) == 0: + sys.exit('no such file: ' + files + '*.tet.dat') + +# get number and name of timesteps +ts = glob.glob(files + '*.tet.dat') + +# loop over timesteps +for time in ts: + + t_step = (time.split('-')[1]).split('.')[0] + file_ts = '%s-%s'%(files,t_step) + + ntelem = 0 + ntpts = 0 + elems = [] + datal = [] + + for file in glob.iglob(file_ts + '*.tet.dat'): + + print 'reading ' + file + fi = open(file, 'r') + + title = fi.readline().split('"')[1].strip() + + vars = fi.readline().split()[2:] + + for i, var in enumerate(vars): + vars[i] = var.strip('"') + + nvar = len(vars) + s = fi.readline().split() + npts = int(s[2]) + nelem = int(s[4]) + + for i in range(npts): + datal.append(fi.readline().split()) + + for i in range(nelem): + elems.append((np.array(fi.readline().split(), dtype=int) - 1 + ntpts).tolist()) + + ntelem += nelem + ntpts += npts + + fi.close() + + + data = np.array(datal, dtype=float).T + elems = np.array(elems, dtype=int).T + + fo = open(file + '.vtu', 'w') + + fo.write('\n') + fo.write('\n') + fo.write('\n') + fo.write('\n') + fo.write('\n') + + for i in range(ntpts): + fo.write(' ' + str('{: .16E}'.format(data[0,i])) + ' ' + str('{: .16E}'.format(data[1,i])) + ' ' + str('{: .16E}'.format(data[2,i])) + '\n') + + fo.write('\n') + fo.write('\n') + fo.write('\n') + fo.write('\n') + + for i in range(ntelem): + fo.write(' ' + str('{:6d}'.format(elems[0,i])) + ' ' + str('{:6d}'.format(elems[1,i])) + ' ' + str('{:6d}'.format(elems[2,i])) + ' ' + str('{:6d}'.format(elems[3,i])) + '\n') + + fo.write('') + + fo.write('\n') + + for i in range(ntelem): + fo.write(' ' + str('{:6d}'.format((i+1)*4)) + '\n') + + fo.write('\n') + + fo.write('\n') + + + for i in range(ntelem): + fo.write(' ' + str('{:6d}'.format(10)) + '\n') + + fo.write('\n') + + fo.write('\n') + + fo.write('\n') + + + for j, var in enumerate(vars[3:]): + fo.write('\n') + + for i in range(npts): + fo.write(' ' + str('{: .16E}'.format(data[j+3,i])) + '\n') + + fo.write('\n') + + + fo.write('\n') + fo.write('\n') + fo.write('\n') + fo.write('\n') + fo.close() + print file + '.vtu written' diff --git a/tools/visz_fine.py b/tools/visz_fine.py new file mode 100644 index 0000000..6db7545 --- /dev/null +++ b/tools/visz_fine.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Simon Kremers (kremers AT geophysik.uni-muenchen.de) +# @author Christian Pelties (pelties AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/pelties) +# @author Tobias Megies (tobias.megies AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/megies) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# +# visualize Seissol output + +import numpy as np +import sys +import matplotlib.pyplot as plt + +try: + infile = sys.argv[1] +except: + raise NameError('no input file!') + +# read file header +f = open(infile) +line1 = f.readline() +line2 = f.readline() +f.close() +variables = line1.replace('"', '').split()[2:] +num_x = int(line2.split()[2]) +num_y = int(line2.split()[4]) +num_var = len(variables) +titl = variables +titl.remove("x") +titl.remove("y") +num_titl = len(titl) + +msg = ["%s: %s" % (i+1, item) for i, item in enumerate(titl)] +msg = "What component do you want to look at?\n" + "\n".join(msg) + "\n" + +comp = -1 +while not 0 <= comp < num_titl: + try: + comp = int(raw_input(msg)) - 1 + except: + pass + +# read in data matrix +data = np.loadtxt(infile, skiprows=2).T +x = data[0].reshape((num_y, num_x)) +y = data[1].reshape((num_y, num_x)) +data = data[2:] +z = data[comp].reshape((num_y, num_x)) +z = np.flipud(z) + +plt.figure() +ax = plt.gca() +ax.set_aspect('equal') +plt.imshow(z, extent=(x.min(), x.max(), y.min(), y.max()), interpolation="nearest") +cb = plt.colorbar() +cb.set_label(titl[comp]) +#plt.triplot(x, y, triangles) +plt.title('%s component of wavefield' % titl[comp]) +plt.xlabel('X [m]') +plt.ylabel('Y [m]') + +plt.show() diff --git a/tools/visz_snap.py b/tools/visz_snap.py new file mode 100644 index 0000000..3c29fb4 --- /dev/null +++ b/tools/visz_snap.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python +## +# @file +# This file is part of SeisSol. +# +# @author Simon Kremers (kremers AT geophysik.uni-muenchen.de) +# @author Christian Pelties (pelties AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/pelties) +# @author Tobias Megies (tobias.megies AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/megies) +# +# @section LICENSE +# Copyright (c) SeisSol Group +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# 3. Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from this +# software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# @section DESCRIPTION +# visualize Seissol output + +import numpy as np +import sys +import matplotlib.pyplot as plt +from StringIO import StringIO + +try: + infile = sys.argv[1] +except: + raise NameError('no input file!') + +# read file header +lines = open(infile).readlines() +line = lines.pop(0) +info = line[line.find('"')+1:line.rfind('"')] + +line = lines.pop(0) +variables = line.replace('"', '').split()[2:] +num_var = len(variables) +titl = variables +titl.remove("x") +titl.remove("y") +num_titl = len(titl) + +line = lines.pop(0) +num_cells = int(line.split()[2]) + +msg = ["%s: %s" % (i+1, item) for i, item in enumerate(titl)] +msg = "What component do you want to look at?\n" + "\n".join(msg) + "\n" + +comp = -1 +while not 0 <= comp < num_titl: + try: + comp = int(raw_input(msg)) - 1 + except: + pass + +# read in data matrix +data = "".join(lines[:num_cells]) +data = StringIO(data) +data = np.loadtxt(data, dtype="float").T + +x = data[0] +y = data[1] +data = data[2:] + +# read in connectivity matrix +triangles = "".join(lines[num_cells:]) +triangles = StringIO(triangles) +triangles = np.loadtxt(triangles, dtype="int") +# matplotlib expects the first point to have index 0, +# in the file it has index 1! +triangles -= 1 + +plt.figure() +ax = plt.gca() +ax.set_aspect('equal') +plt.tripcolor(x, y, triangles, data[comp, :], shading='faceted') +cb = plt.colorbar() +cb.set_label(titl[comp]) +#plt.triplot(x, y, triangles) +plt.title('%s -- %s component of wavefield' % (info, titl[comp])) +plt.xlabel('X [m]') +plt.ylabel('Y [m]') + +plt.show()