2016-PippingKornhuberRosenauOncken to 75211655b5f958239e3a427827f0a4f56512eb28 · podlesny / dune-tectonic

Some changes are not shown.

For a faster browsing experience, only 110 of 357 files are shown. Download one of the files below to see all changes.

CMakeLists.txt

+1 −1

Original line number	Diff line number	Diff line
		@@ -17,7 +17,7 @@ include(DuneMacros)

		# start a dune project with information from dune.module
		dune_project()

		dune_enable_all_packages()
		find_package(HDF5 COMPONENTS C REQUIRED)

		add_subdirectory("src")

data/2d-dip-contours-performance.tikz

0 → 100644

+91 −0

Original line number	Diff line number	Diff line
		\pgfplotsarraynew\heights{%
		-40\\-20\\-10\\-5\\-2.5\\%
		0\\%
		2.5\\5\\10\\20\\40\\%
		}
		\def\fnamebase{2d-dip-contours}
		\def\fname{rfpitol=100e-7}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:threequakes:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:threequakes:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		xmin=\threequakesmin, xmax=\threequakesmax, max space between ticks=40pt,
		x tick label style={ /pgf/number format/1000 sep={} },
		tick label style={font=\footnotesize},
		label style={font=\small},
		group style={
		x descriptions at = edge bottom,
		group size=1 by 4,
		vertical sep=0cm
		},
		width=12cm,
		enlargelimits=false]
		\nextgroupplot[
		semithick,
		height = 6.5cm,
		/pgf/number format/1000 sep={},
		ymax=0.7,
		colormap/jet,
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small,
		at={(1.05,1)},
		anchor=north west,
		fill=none},
		ylabel = distance from trench, y unit = m,
		legend cell align=right,
		contour/labels=false,
		contour prepared]
		\pgfplotsinvokeforeach{1,...,9} { % level 0 and 10 are empty
		\pgfplotscolormapaccess[0:10]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\def\fnameX{generated/\fnamebase:\fname:level:%
		\pgfplotsarrayvalueofelem#1\of\heights.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\heights}{\micro\meter}}
		};
		\nextgroupplot[
		semithick,
		height = 3.5cm,
		ylabel style={align=center}, ylabel = time step\\size, y unit = s,
		ytick={1e-3,1e-2,1e-1},
		ymax = 1, ymin = 1e-4,
		ymode = log]
		\addplot[const plot, mark=none] table[col sep=comma, y=timeIncrement]
		{generated/2d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		height = 3.5cm,
		ylabel style={align=center}, ylabel=fixed-point\\iterations,
		ytick={2,4,6,8},
		ymin=0, ymax=10]
		\addplot[const plot, mark=none] table[col sep=comma, y=fixedPointIterations]
		{generated/2d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		height = 3.5cm,
		xlabel = time, x unit = s,
		ylabel style={align=center}, ylabel=multigrid\\iterations,
		ytick={5,10,15,20},
		ymin=0, ymax=25]
		\addplot[const plot, mark=none] table[col sep=comma, y=multiGridIterations]
		{generated/2d-performance:\fname.csv};
		\end{groupplot}
		\end{tikzpicture}

data/2d-dip-single-points.tikz

0 → 100644

+51 −0

Original line number	Diff line number	Diff line
		\def\fnamebase{dip-single-points}
		\def\fname{rfpitol=100e-7}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:threequakes:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:threequakes:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		xmin=\threequakesmin, xmax=\threequakesmax, max space between ticks=40pt,
		/pgf/number format/1000 sep={},
		x tick label style={
		/pgf/number format/.cd,
		1000 sep={},
		/tikz/.cd
		},
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small, at={(1.05,1)},
		anchor=north west,
		fill=none},
		legend entries={
		\SI{15}{\centi\meter},
		\SI{30}{\centi\meter},
		\SI{45}{\centi\meter}
		},
		ylabel style={align=center}, ylabel=vertical surface\\displacement, y unit = m,
		change y base, y SI prefix=micro,
		xlabel = time, x unit=s,
		width = 12cm,
		height = 4cm,
		semithick]
		\addplot[width=2pt] table[col sep=comma, x index = 0, y index=1]
		{generated/\fnamebase:\fname.csv};
		\addplot[width=2pt, dashed] table[col sep=comma, x index = 0, y index=2]
		{generated/\fnamebase:\fname.csv};
		\addplot[width=2pt, dotted] table[col sep=comma, x index = 0, y index=3]
		{generated/\fnamebase:\fname.csv};
		\end{axis}
		\end{tikzpicture}

data/2d-effort-over-tolerance.tikz

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		semithick,
		width = 12cm,
		height = 4cm,
		xmode=log,
		%ymin = 0,
		% scaled y ticks=base 10:0,
		xlabel = fixed point tolerance,
		ylabel style={align=center}, ylabel = multigrid\\iterations, % (multigrid steps)
		extra x ticks = 1e-5,
		extra x tick labels = { time-stepping\\tolerance },
		extra x tick style = { align=center, grid = major, ticklabel pos=right }
		]
		\addplot[mark=+] table[col sep=comma, y=mg] {generated/fpi-data.csv};
		\end{axis}
		\end{tikzpicture}

data/2d-velocity-contours-threequakes.tikz

0 → 100644

+61 −0

Original line number	Diff line number	Diff line
		\pgfplotsarraynew\contourlevels{%
		1\\3\\
		10\\30\\
		100\\300\\
		1000\\3000\\
		}

		\def\fname{threequakes:rfpitol=100e-7}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		xmin=\threequakesmin, xmax=\threequakesmax, max space between ticks=40pt,
		ymax=0.85,
		colormap/jet,
		x tick label style={
		/pgf/number format/.cd,
		1000 sep={},
		/tikz/.cd
		},
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small,
		at={(1.05,1)},
		anchor=north west,
		fill=none},
		ylabel = distance from trench, y unit = m,
		xlabel = time, x unit=s,
		width = 12cm, height = 6.5cm,
		legend cell align=right,
		contour/labels=false,
		enlargelimits=false,
		semithick]

		\pgfplotsinvokeforeach{0,2,4,6} {
		\pgfplotscolormapaccess[0:7]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\def\fnameX{generated/2d-velocity-contours:\fname:level:%
		\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\contourlevels}{\micro\meter/\second}}
		};
		\end{axis}
		\end{tikzpicture}

data/2d-velocity-contours-zoom.tikz

0 → 100644

+51 −0

Original line number	Diff line number	Diff line
		\pgfplotsarraynew\contourlevels{%
		1\\3\\
		10\\30\\
		100\\300\\
		1000\\3000\\
		}

		\def\fname{zoom:rfpitol=100e-7}
		\begin{tikzpicture}[trim axis left, trim axis right]
		\begin{axis}[
		max space between ticks=40pt,
		ymax=0.7,
		colormap/jet,
		x tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		1000 sep={},
		/tikz/.cd
		},
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		legend style={font=\small,
		at={(1.05,1)},
		anchor=north west,
		fill=none},
		ylabel = distance from trench, y unit = m,
		xlabel = time, x unit=s,
		width = 12cm, height = 6.5cm,
		legend cell align=right,
		enlargelimits=false,
		contour/labels=false,
		semithick]

		\pgfplotsinvokeforeach{0,...,6} { % level 7 is empty
		\pgfplotscolormapaccess[0:7]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\def\fnameX{generated/2d-velocity-contours:\fname:level:%
		\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\contourlevels}{\micro\meter/\second}}
		};
		\end{axis}
		\end{tikzpicture}

data/3d-performance.tikz

0 → 100644

+57 −0

Original line number	Diff line number	Diff line
		\def\fname{rtol=1e-5_diam=1e-2}

		\newread\threequakesminread
		\openin\threequakesminread=generated/timeframe:min:threequakes:\fname.tex
		\read\threequakesminread to \threequakesmin
		\closein\threequakesminread

		\newread\threequakesmaxread
		\openin\threequakesmaxread=generated/timeframe:max:threequakes:\fname.tex
		\read\threequakesmaxread to \threequakesmax
		\closein\threequakesmaxread

		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		xmin=\threequakesmin, xmax=\threequakesmax,
		max space between ticks=40pt,
		x tick label style={
		/pgf/number format/.cd,
		1000 sep={},
		/tikz/.cd
		},
		tick label style={font=\footnotesize},
		label style={font=\small},
		group style={
		x descriptions at = edge bottom,
		group size=1 by 3,
		vertical sep=0cm
		},
		height=3.5cm, width=12cm,
		enlargelimits=false]

		\nextgroupplot[
		semithick,
		ylabel style={align=center}, ylabel = time step\\size, y unit = s,
		ytick={1e-3,1e-2,1e-1},
		ymax = 1, ymin = 1e-4,
		ymode = log]
		\addplot[const plot, mark=none] table[col sep=comma, y=timeIncrement]
		{generated/3d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		ylabel style={align=center}, ylabel = fixed-point\\iterations,
		ytick={2,4,6,8},
		ymin=0, ymax=10]
		\addplot[const plot, mark=none] table[col sep=comma, y=fixedPointIterations]
		{generated/3d-performance:\fname.csv};
		\nextgroupplot[
		semithick,
		xlabel = time, x unit = s,
		ylabel style={align=center}, ylabel=multigrid\\iterations,
		ytick={5,10,15,20},
		ymin=0, ymax=25
		]
		\addplot[const plot, mark=none] table[col sep=comma, y=multiGridIterations]
		{generated/3d-performance:\fname.csv};
		\end{groupplot}
		\end{tikzpicture}

data/3d-velocity-contours.tikz

0 → 100644

+156 −0

Original line number	Diff line number	Diff line
		\def\fname{generated/3d-velocity-contours:rtol=1e-5_diam=1e-2}
		\pgfplotstableread[col sep=comma]{\fname:times.csv}\myloadedtable

		\pgfplotsarraynew\contourlevels{%
		1\\2\\3\\5\\%
		10\\20\\30\\50\\%
		100\\200\\300\\500\\%
		1000\\
		}

		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		group style={
		y descriptions at = edge left,
		group size=6 by 1,
		horizontal sep=0cm
		},
		ymin=-0.30, ymax= 0.30,
		enlarge x limits=false,
		colormap/jet,
		y tick label style={
		/pgf/number format/.cd,
		fixed, fixed zerofill, precision=2,
		/tikz/.cd
		},
		ytick = {-0.30,-0.20,-0.10,0.00,0.10,0.20,0.30},
		tick label style={font=\footnotesize},
		label style={font=\small},
		width = 3.5cm, height = 6cm,
		enlargelimits=false,
		contour/labels=false,
		%
		groupplot xlabel={distance from trench [\si{\meter}]},
		]

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{0}{times}\of\myloadedtable%
		$t_0 \approx
		\SI[round-mode=places,round-precision=0]{\pgfplotsretval}{\second}$
		},
		ylabel = width, y unit = m]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} {
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}

		\pgfplotstablegetelem{0}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{1}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{1}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{2}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{2}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{3}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		},
		legend columns=4,
		legend cell align=right,
		legend style={font=\small,
		at={(0,1.05)},
		anchor=south,
		fill=none},
		]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{3}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}, forget plot] table \fnameX;
		\addlegendimage{line legend,color=mycolor#1}
		\addlegendentry{%
		\SI{\pgfplotsarrayvalueofelem#1\of\contourlevels}{\micro\meter/\second}}
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{4}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{4}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};

		\nextgroupplot[semithick, xlabel = {
		\pgfplotstablegetelem{5}{timeOffsets}\of\myloadedtable%
		$t_0 +
		\SI[round-mode=places,round-precision=2]{\pgfplotsretval}{\second}$
		}]
		\draw[color=gray!20] (0.162533,-0.30) -- (0.162533,+0.30); % X
		\draw[color=gray!20] (0.362533+0.05,-0.30) -- (0.362533-0.05,+0.30); % Y
		\pgfplotsinvokeforeach{0,...,12} { % level 13 and 14 are empty
		\pgfplotscolormapaccess[0:14]{#1}{jet}
		\def\TEMP{\definecolor{mycolor#1}{rgb}}
		\expandafter\TEMP\expandafter{\pgfmathresult}
		\pgfplotstablegetelem{5}{timeSteps}\of\myloadedtable
		\edef\fnameX{\fname:%
		step:\pgfplotsretval:%
		level:\pgfplotsarrayvalueofelem#1\of\contourlevels.tex}
		\addplot[contour prepared={draw color=mycolor#1}] table \fnameX;
		};
		\end{groupplot}
		\end{tikzpicture}

data/boxplot.tikz

0 → 100644

+43 −0

Original line number	Diff line number	Diff line
		\def\simulationtag{rfpitol=100e-7}
		\begin{tikzpicture}[trim axis group left, trim axis group right]
		\begin{groupplot}[
		tick label style={font=\footnotesize},
		label style={font=\small},
		%
		group style={
		y descriptions at = edge left,
		group size=3 by 1,
		horizontal sep=0.75cm
		},
		height=4cm,
		width=4.5cm,
		%
		ytick={1,2,3},
		yticklabels={experiment, simulation}
		]
		%
		\nextgroupplot[semithick, xlabel = recurrence time, x unit = s, xmode=log,
		log ticks with fixed point, xtick={5,10,20,40}]
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:lab:recurrence.tex};
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:simulation:\simulationtag:recurrence.tex};
		%
		\nextgroupplot[semithick, xlabel = rupture width, x unit=m, xmin=0, xmax=0.4,
		extra x ticks = 0.2,
		extra x tick labels = ,
		extra x tick style = { grid = major }
		]
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:lab:ruptureWidth.tex};
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:simulation:\simulationtag:ruptureWidth.tex};
		%
		\nextgroupplot[semithick, xlabel = peak slip, x unit=mm, xmode=log, log ticks with fixed point,
		xtick={0.03, 0.06, 0.12}]
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:lab:peakSlip.tex};
		\addplot[mark=+, boxplot={ box extend = 0.5 }]
		table[y index=0] {generated/boxplot-data:simulation:\simulationtag:peakSlip.tex};
		\end{groupplot}
		\end{tikzpicture}

data/includes.tex

0 → 100644

+51 −0

Original line number	Diff line number	Diff line
		\usepackage{pgfplots}
		\pgfplotsset{compat=1.11} % FIXME: 1.12 would be nice; debian:8.7 only has 1.11
		\usepackage{pgfplotstable}
		\usepgfplotslibrary{groupplots}
		\usepgfplotslibrary{statistics}
		\usepgfplotslibrary{units}


		%% Typeset the mu from '[xy] SI prefix=micro' as an upright mu
		%% From https://tex.stackexchange.com/a/224574
		\pgfplotsset{
		x SI prefix/micro/.style={/pgfplots/axis base prefix={axis x base 6 prefix \micro}},
		y SI prefix/micro/.style={/pgfplots/axis base prefix={axis y base 6 prefix \micro}},
		z SI prefix/micro/.style={/pgfplots/axis base prefix={axis z base 6 prefix \micro}},
		unit code/.code 2 args={\si{#1#2}},
		}

		%% Add support for 'groupplot [xy]label'
		%% From http://tex.stackexchange.com/a/117935/16940, see also
		%% https://sourceforge.net/p/pgfplots/feature-requests/48/
		\makeatletter
		\pgfplotsset{
		groupplot xlabel/.initial={},
		every groupplot x label/.style={
		at={($({\pgfplots@group@name\space c1r\pgfplots@group@rows.west}\|-{\pgfplots@group@name\space c1r\pgfplots@group@rows.outer south})!0.5!({\pgfplots@group@name\space c\pgfplots@group@columns r\pgfplots@group@rows.east}\|-{\pgfplots@group@name\space c\pgfplots@group@columns r\pgfplots@group@rows.outer south})$)},
		anchor=north,
		},
		groupplot ylabel/.initial={},
		every groupplot y label/.style={
		rotate=90,
		at={($({\pgfplots@group@name\space c1r1.north}-\|{\pgfplots@group@name\space c1r1.outer
		west})!0.5!({\pgfplots@group@name\space c1r\pgfplots@group@rows.south}-\|{\pgfplots@group@name\space c1r\pgfplots@group@rows.outer west})$)},
		anchor=south
		},
		execute at end groupplot/.code={%
		\node [/pgfplots/every groupplot x label]
		{\pgfkeysvalueof{/pgfplots/groupplot xlabel}};
		\node [/pgfplots/every groupplot y label]
		{\pgfkeysvalueof{/pgfplots/groupplot ylabel}};
		}
		}

		\def\endpgfplots@environment@groupplot{%
		\endpgfplots@environment@opt%
		\pgfkeys{/pgfplots/execute at end groupplot}%
		\endgroup%
		}
		\makeatother

		%% Have \includegraphics{} support tikz files
		\usepackage{tikzscale}
		No newline at end of file

data/standalone/2d-dip-contours-performance.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-dip-contours-performance.tikz}
		\end{document}

data/standalone/2d-dip-single-points.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-dip-single-points.tikz}
		\end{document}

data/standalone/2d-effort-over-tolerance.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-effort-over-tolerance.tikz}
		\end{document}

data/standalone/2d-velocity-contours-threequakes.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-velocity-contours-threequakes}
		\end{document}

data/standalone/2d-velocity-contours-zoom.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{2d-velocity-contours-zoom}
		\end{document}

data/standalone/3d-performance.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{3d-performance.tikz}
		\end{document}

data/standalone/3d-velocity-contours.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\includegraphics{3d-velocity-contours.tikz}
		\end{document}

data/standalone/Makefile

0 → 100644

+18 −0

Original line number	Diff line number	Diff line
		PDFs=\
		2d-dip-contours-performance.pdf \
		2d-dip-single-points.pdf \
		2d-effort-over-tolerance.pdf \
		2d-velocity-contours-threequakes.pdf \
		2d-velocity-contours-zoom.pdf \
		3d-performance.pdf \
		3d-velocity-contours.pdf \
		boxplot.pdf
		PNGs=$(PDFs:.pdf=.png)

		pdf: $(PDFs)
		%.pdf: standalone/%.tex
		latexmk -pdf $<

		png: $(PNGs)
		%.png: %.pdf
		convert -density 300 $< -quality 100 $@

data/standalone/boxplot.tex

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		\documentclass[border={3cm 0cm}]{standalone}
		\usepackage{siunitx}
		\input{includes}

		\begin{document}
		\input{boxplot.tikz}
		\end{document}

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

0 → 100644

+108 −0

Original line number	Diff line number	Diff line
		source('tools/support/findQuakes.R')
		source('tools/support/writeContours.R')
		Rcpp::sourceCpp('tools/support/trapezoidal.cpp')

		finalTime <- 1000 # s
		specialTrenchDistances <- c(0.15,0.30,0.45) # m
		convergenceVelocity <- 5e-5 # m/s

		last <- function(x) x[[length(x)]]
		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (basedir in c("rfpitol=100e-7")) {
		dir <- file.path(directories[['simulation']],
		'2d-lab-fpi-tolerance', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
		relativeTime <- h5file['relativeTime'][]
		realTime <- finalTime * relativeTime

		calcMask <- relativeTime >= 0.7
		calcTime <- realTime[calcMask]
		calcDuration <- last(calcTime) - calcTime[[1]]
		calcRange <- range(which(calcMask))
		calcLength <- sum(calcMask)

		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a very low velocity here.
		quakes <- findQuakes(1e-6 + convergenceVelocity,
		h5file['/frictionalBoundary/velocity'][,,],
		indices = calcRange[1]:calcRange[2])
		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]
		quakes$duration <- quakes$ending - quakes$beginning
		numQuakes <- nrow(quakes)

		relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
		quakes[[numQuakes, 'ending']]), f=0.02)
		threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])

		plotMask <- threeQuakeTimeMask
		plotTime <- realTime[plotMask]
		plotRange <- range(which(plotMask))
		plotLength <- sum(plotMask)

		write(relaxedTime[[1]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'min', 'threequakes', basedir),
		'tex')))
		write(relaxedTime[[2]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'max', 'threequakes', basedir),
		'tex')))

		surfaceCoordinates <- h5file['/surface/coordinates'][]
		surfaceTrenchDistance <- surfaceCoordinates[,1] / cos(atan(.27))
		perm <- order(surfaceTrenchDistance)

		displacement <- h5file['/surface/displacement']
		# This is the displacement imposed through the Dirichlet condition at the last node
		displacementOffset <- displacement[calcRange[1]:calcRange[2],last(perm),1]
		balancedSurfaceDisplacement <- matrix(nrow = plotLength, ncol = displacement@dim[2])
		for (k in 1:displacement@dim[2]) {
		d <- displacement[calcRange[1]:calcRange[2],k,1:2]
		d[,,1] <- d[,,1] - displacementOffset
		# We're in a tilted coordinate system
		dv <- -sin(atan(.27))d[,,1] + cos(atan(.27))d[,,2]
		meanVertSurfaceDisplacement <- trapezoidal(calcTime, dv) / calcDuration
		balancedSurfaceDisplacement[,k] <- (dv - meanVertSurfaceDisplacement)[
		(plotRange[1] - calcRange[1] + 1):(plotRange[2] - calcRange[1] + 1)
		]
		}

		## Interpolate velocity to special points from surrounding velocities
		{
		data <- matrix(nrow=plotLength, ncol=1+length(specialTrenchDistances))
		data[,1] <- plotTime
		for (k in seq(plotLength)) {
		interpolant <- approxfun(surfaceTrenchDistance[perm],
		balancedSurfaceDisplacement[k,perm])
		for (i in seq(specialTrenchDistances)) {
		specialTrenchDistance <- specialTrenchDistances[[i]]
		data[k,i+1] <- interpolant(specialTrenchDistance)
		}
		}
		singleDataName <- file.path(directories[['output']],
		paste.(pasteColon('dip-single-points', basedir),
		'csv'))
		write.csv(data, singleDataName, row.names = FALSE)
		}
		h5::h5close(h5file)

		printlevels <- c('-40','-20','-10','-5','-2.5','0','2.5','5','10','20','40')
		levels <- as.numeric(printlevels) * 1e-6

		ret <- contourLines(plotTime,
		surfaceTrenchDistance[perm],
		balancedSurfaceDisplacement[,perm],
		levels=levels)

		for (i in seq(printlevels))
		writeContours(ret, levels[[i]],
		file.path(directories[['output']],
		paste.(pasteColon('2d-dip-contours', basedir,
		'level', printlevels[[i]]),
		'tex')))
		}

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

0 → 100644

+68 −0

Original line number	Diff line number	Diff line
		source('tools/support/findQuakes.R')

		finalTime <- 1000 # s
		convergenceVelocity <- 5e-5 # m/s
		criticalVelocity <- 1000e-6 + convergenceVelocity

		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (basedir in c("rfpitol=100e-7")) {
		dir <- file.path(directories[['simulation']],
		'2d-lab-fpi-tolerance', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
		relativeTime <- h5file['relativeTime'][]
		realTime <- finalTime * relativeTime

		velocityProxy<- h5file['/frictionalBoundary/velocity']

		quakes <- findQuakes(criticalVelocity, velocityProxy,
		indices = 1:dim(velocityProxy)[1], 1)

		basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]

		maximumVelocities<- maxVelocity(velocityProxy[,,1])
		displacement <- h5file['/frictionalBoundary/displacement']
		numVertices <- dim(displacement)[[2]]
		for (quakeID in seq(nrow(quakes))) {
		qb <- quakes[[quakeID,'beginningIndex']]
		qe <- quakes[[quakeID,'endingIndex']]
		localSlippingTimes <- abs(velocityProxy[qb:qe,,1][,,1]) > criticalVelocity

		qd <- displacement[qb:qe,,1]
		slip <- vector(mode='numeric', length=numVertices)
		for (i in seq(numVertices)) {
		if (any(localSlippingTimes[,i])) {
		begs <- positiveStarts(localSlippingTimes[,i])
		ends <- negativeStarts(localSlippingTimes[,i])
		slip[i]<- sum(qd[ends,i,] - qd[begs,i,])
		}
		}
		quakes[quakeID,'peakSlip'] <- max(abs(slip))
		quakes[quakeID,'peakSlipRate'] <- max(maximumVelocities[qb:qe])

		maxRuptureWidth <- 0
		for (ts in seq(dim(localSlippingTimes)[[1]])) {
		st <- localSlippingTimes[ts,]
		if (!any(st))
		next;

		slippingCoordinates <- basalCoordinates[st,1] # x-coordinate only
		maxRuptureWidth <- max(maxRuptureWidth, diff(range(slippingCoordinates)))
		}
		quakes[quakeID,'ruptureWidth'] <- maxRuptureWidth
		}
		h5::h5close(h5file)

		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]

		write.csv(quakes[c('beginning','ending',
		'peakSlip','peakSlipRate',
		'ruptureWidth')],
		row.names = FALSE, quote = FALSE,
		file = file.path(directories[['output']],
		paste.(pasteColon('events', basedir), 'csv')))
		}

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

0 → 100644

+41 −0

Original line number	Diff line number	Diff line
		rfpitols <- c('1e-7', '2e-7', '3e-7', '5e-7',
		'10e-7', '20e-7', '30e-7', '50e-7',
		'100e-7', '200e-7', '300e-7', '500e-7',
		'1000e-7', '2000e-7', '3000e-7', '5000e-7',
		'10000e-7', '20000e-7', '30000e-7', '50000e-7',
		'100000e-7')
		numEntries <- length(rfpitols)

		data <- data.frame(row.names = rfpitols,
		tol = rep(NA, numEntries),
		time = rep(NA, numEntries),
		fpi = rep(NA, numEntries),
		mg = rep(NA, numEntries))

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (rfpitol in rfpitols) {
		basedir <- paste('rfpitol', rfpitol, sep='=')
		dir <- file.path(directories[['simulation']],
		'2d-lab-fpi-tolerance', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')

		data[rfpitol,'tol'] <- as.numeric(rfpitol)

		relativeTimeProxy <- h5file['/relativeTime']
		relativeTimeLen <- relativeTimeProxy@dim
		data[rfpitol,'time'] <- relativeTimeProxy[relativeTimeProxy@dim]

		## FIXME: why do we drop the first entry?
		fixedPointIterationsProxy <- h5file["/iterations/fixedPoint/total"]
		fixedPointIterationsLen <- fixedPointIterationsProxy@dim
		data[rfpitol,'fpi'] <- sum(fixedPointIterationsProxy[2:fixedPointIterationsLen])

		multiGridIterationsProxy <- h5file["/iterations/multiGrid/total"]
		multiGridIterationsLen <- multiGridIterationsProxy@dim
		data[rfpitol,'mg'] <- sum(multiGridIterationsProxy[2:multiGridIterationsLen])
		h5::h5close(h5file)
		}

		write.csv(data, file.path(directories[['output']], 'fpi-data.csv'),
		row.names = FALSE, quote = FALSE)

data/tools/2d-performance.R

0 → 100644

+58 −0

Original line number	Diff line number	Diff line
		source('tools/support/findQuakes.R')

		finalTime <- 1000 # s
		convergenceVelocity <- 5e-5 # m/s

		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (basedir in c("rfpitol=100e-7")) {
		dir <- file.path(directories[['simulation']],
		'2d-lab-fpi-tolerance', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
		relativeTime <- h5file['relativeTime'][]
		realTime <- finalTime * relativeTime

		velocityProxy<- h5file['/frictionalBoundary/velocity']

		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a very low velocity here.
		quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
		indices = 1:dim(velocityProxy)[1], 1)
		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]
		quakes$duration <- quakes$ending - quakes$beginning
		numQuakes <- nrow(quakes)

		relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
		quakes[[numQuakes, 'ending']]), f=0.02)
		threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
		plotMask <- threeQuakeTimeMask
		plotIndices<- which(plotMask)
		plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)

		write(relaxedTime[[1]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'min', 'threequakes',
		basedir), 'tex')))
		write(relaxedTime[[2]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'max', 'threequakes',
		basedir), 'tex')))

		timeWindow = realTime[plotIndices]

		relativeTau <- h5file['relativeTimeIncrement'][]
		fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
		multiGridIterations <- h5file['/iterations/multiGrid/final'][]
		write.csv(data.frame(time = timeWindow,
		timeIncrement = finalTime * relativeTau[plotIndices],
		fixedPointIterations = fixedPointIterations[plotIndices],
		multiGridIterations = multiGridIterations[plotIndices]),
		file.path(directories[['output']],
		paste.(pasteColon('2d-performance', basedir), 'csv')),
		row.names = FALSE, quote = FALSE)
		h5::h5close(h5file)
		}

data/tools/2d-performance.py

0 → 100644

+76 −0

Original line number	Diff line number	Diff line
		import ConfigParser as cp
		import os
		import numpy as np
		import csv
		import h5py

		from support.find_quakes import find_quakes

		NBODIES = 2
		FINAL_TIME = 1000 # s
		FINAL_VELOCITY = 1e-5 # m/s
		THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY

		TANGENTIAL_COORDS = 1

		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('config.ini')
		config.read(config_path)
		sim_path = config.get('directories', 'simulation')
		exp_path = config.get('directories', 'experiment')
		out_path = config.get('directories', 'output')

		# create output directory
		out_dir = os.path.join(out_path)
		if not os.path.exists(out_dir):
		os.mkdir(out_dir)

		h5path = os.path.join(sim_path)
		h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		# read time
		relative_time = np.array(h5file['relativeTime'])
		real_time = relative_time * FINAL_TIME


		velocityProxy<- h5file['/frictionalBoundary/velocity']

		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a very low velocity here.
		quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
		indices = 1:dim(velocityProxy)[1], 1)
		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]
		quakes$duration <- quakes$ending - quakes$beginning
		numQuakes <- nrow(quakes)

		relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
		quakes[[numQuakes, 'ending']]), f=0.02)
		threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
		plotMask <- threeQuakeTimeMask
		plotIndices<- which(plotMask)
		plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)

		write(relaxedTime[[1]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'min', 'threequakes',
		basedir), 'tex')))
		write(relaxedTime[[2]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'max', 'threequakes',
		basedir), 'tex')))

		timeWindow = realTime[plotIndices]

		relativeTau <- h5file['relativeTimeIncrement'][]
		fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
		multiGridIterations <- h5file['/iterations/multiGrid/final'][]
		write.csv(data.frame(time = timeWindow,
		timeIncrement = finalTime * relativeTau[plotIndices],
		fixedPointIterations = fixedPointIterations[plotIndices],
		multiGridIterations = multiGridIterations[plotIndices]),
		file.path(directories[['output']],
		paste.(pasteColon('2d-performance', basedir), 'csv')),
		row.names = FALSE, quote = FALSE)
		h5::h5close(h5file)

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

0 → 100644

+95 −0

Original line number	Diff line number	Diff line
		source('tools/support/findQuakes.R')
		source('tools/support/writeContours.R')

		finalTime <- 1000 # s
		convergenceVelocity <- 5e-5 # m/s

		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (basedir in c("rfpitol=100e-7")) {
		dir <- file.path(directories[['simulation']],
		'2d-lab-fpi-tolerance', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
		relativeTime <- h5file['relativeTime'][]
		realTime <- finalTime * relativeTime

		velocityProxy<- h5file['/frictionalBoundary/velocity']

		basalCoordinates <- h5file['/frictionalBoundary/coordinates'][]
		basalTrenchDistance <- basalCoordinates[,1]
		perm <- order(basalTrenchDistance)
		sortedBasalTrenchDistance <- basalTrenchDistance[perm]

		{
		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a very low velocity here.
		quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
		indices = 1:dim(velocityProxy)[1], 1)
		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]
		quakes$duration <- quakes$ending - quakes$beginning
		numQuakes <- nrow(quakes)

		relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
		quakes[[numQuakes, 'ending']]), f=0.02)
		plotMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
		plotIndices<- which(plotMask)

		write(relaxedTime[[1]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'min', 'threequakes',
		basedir), 'tex')))
		write(relaxedTime[[2]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'max', 'threequakes',
		basedir), 'tex')))

		printlevels <- c('1000','100','10','1')
		levels <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
		ret <- contourLines(realTime[plotIndices],
		sortedBasalTrenchDistance,
		abs(velocityProxy[plotIndices,perm,1][,,1]),
		levels = levels)

		for (i in seq(printlevels))
		writeContours(ret, levels[[i]],
		file.path(directories[['output']],
		paste.(pasteColon('2d-velocity-contours',
		'threequakes', basedir, 'level',
		printlevels[[i]]), 'tex')))
		}
		{
		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a rather low velocity here.
		quakes <- findQuakes(300e-6 + convergenceVelocity, velocityProxy,
		indices = 1:dim(velocityProxy)[1], 1)
		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]
		quakes$duration <- quakes$ending - quakes$beginning
		numQuakes <- nrow(quakes)
		quake <- quakes[numQuakes,]
		relaxedTime <-
		c(quake[['beginning']] - 0.9*(quake[['ending']] - quake[['beginning']]),
		quake[['ending']] + 0.1*(quake[['ending']] - quake[['beginning']]))
		plotMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
		plotIndices<- which(plotMask)

		printlevels <- c('3000','1000','300','100','30','10','3','1')
		levels <- 1e-6 * as.numeric(printlevels) + convergenceVelocity
		ret <- contourLines(realTime[plotIndices],
		sortedBasalTrenchDistance,
		abs(velocityProxy[plotIndices,perm,1][,,1]),
		levels = levels)

		for (i in seq(printlevels))
		writeContours(ret, levels[[i]],
		file.path(directories[['output']],
		paste.(pasteColon('2d-velocity-contours',
		'zoom', basedir, 'level',
		printlevels[[i]]), 'tex')))
		}
		h5::h5close(h5file)
		}

data/tools/2d_event_writer.py

0 → 100644

+112 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import numpy as np
		import csv
		import h5py

		from support.maximum import maximum
		from support.norm import norm
		from support.find_quakes import find_quakes
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings
		from support.max_distance import max_distance


		NBODIES = 2
		FINAL_TIME = 1000 # s
		FINAL_VELOCITY = 1e-5 # m/s
		THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY

		TANGENTIAL_COORDS = 1

		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('config.ini')
		config.read(config_path)
		sim_path = config.get('directories', 'simulation')
		exp_path = config.get('directories', 'experiment')
		out_path = config.get('directories', 'output')

		# create output directory
		out_dir = os.path.join(out_path)
		if not os.path.exists(out_dir):
		os.mkdir(out_dir)

		h5path = os.path.join(sim_path)
		h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		# read time
		relative_time = np.array(h5file['relativeTime'])
		real_time = relative_time * FINAL_TIME

		for body_ID in range(NBODIES):
		body = 'body' + str(body_ID)

		if body not in h5file:
		continue

		# read data
		coordinates = np.array(h5file[body + '/coordinates'])
		displacement = np.array(h5file[body + '/displacement'])
		velocity = np.array(h5file[body + '/velocity'])

		num_vertices = displacement.shape[1]
		velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
		maximum_velocity = maximum(velocity_norm)

		[quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
		maximum_velocity)

		quakes = {}
		for quake_ID in range(len(quake_starts)):
		quake_start = int(quake_starts[quake_ID])
		quake_end = int(quake_ends[quake_ID])

		quake = {}
		quake['start'] = real_time[quake_start]
		quake['end'] = real_time[quake_end]

		local_slipping_times = velocity_norm[quake_start:quake_end, :] \
		> THRESHOLD_VELOCITY
		quake_displacement = displacement[quake_start:quake_end, :, :]
		slip = np.zeros(num_vertices)

		for i in range(num_vertices):
		if any(local_slipping_times[:, i]):
		starts = slip_beginnings(local_slipping_times[:, i])
		ends = slip_endings(local_slipping_times[:, i])
		slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
		- quake_displacement[starts, i, :])

		quake['peakSlip'] = max(slip)
		quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])

		max_rupture_width = 0
		for time_step in range(local_slipping_times.shape[0]):
		slipping_time = local_slipping_times[time_step, :]
		if not any(slipping_time):
		continue

		slipping_coords = coordinates[slipping_time] \
		+ displacement[quake_start + time_step, slipping_time]
		if len(slipping_coords) > 1:
		pair = np.array(max_distance(slipping_coords))
		max_rupture_width = max(max_rupture_width,
		np.linalg.norm(pair[0] - pair[1]))

		quake['ruptureWidth'] = max_rupture_width
		quakes[quake_ID] = quake

		# output quake data to csv file
		csv_columns = quakes[0].keys()
		csv_file_path = os.path.join(out_path, 'events' + str(body_ID) + '.csv')
		try:
		with open(csv_file_path, 'w') as csv_file:
		writer = csv.DictWriter(csv_file, fieldnames=csv_columns)
		writer.writeheader()
		for quake_ID in quakes:
		writer.writerow(quakes[quake_ID])
		except IOError:
		print('I/O error')

		h5file.close()

data/tools/2d_event_writer.py~

0 → 100644

+112 −0

Original line number	Diff line number	Diff line
		import ConfigParser as cp
		import os
		import numpy as np
		import csv
		import h5py

		from support.maximum import maximum
		from support.norm import norm
		from support.find_quakes import find_quakes
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings
		from support.max_distance import max_distance


		NBODIES = 2
		FINAL_TIME = 1000 # s
		FINAL_VELOCITY = 1e-5 # m/s
		THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY

		TANGENTIAL_COORDS = 1

		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('config.ini')
		config.read(config_path)
		sim_path = config.get('directories', 'simulation')
		exp_path = config.get('directories', 'experiment')
		out_path = config.get('directories', 'output')

		# create output directory
		out_dir = os.path.join(out_path)
		if not os.path.exists(out_dir):
		os.mkdir(out_dir)

		h5path = os.path.join(sim_path)
		h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		# read time
		relative_time = np.array(h5file['relativeTime'])
		real_time = relative_time * FINAL_TIME

		for body_ID in range(NBODIES):
		body = 'body' + str(body_ID)

		if body not in h5file:
		continue

		# read data
		coordinates = np.array(h5file[body + '/coordinates'])
		displacement = np.array(h5file[body + '/displacement'])
		velocity = np.array(h5file[body + '/velocity'])

		num_vertices = displacement.shape[1]
		velocity_norm = norm(velocity[:, :, TANGENTIAL_COORDS:])
		maximum_velocity = maximum(velocity_norm)

		[quake_starts, quake_ends] = find_quakes(THRESHOLD_VELOCITY,
		maximum_velocity)

		quakes = {}
		for quake_ID in range(len(quake_starts)):
		quake_start = int(quake_starts[quake_ID])
		quake_end = int(quake_ends[quake_ID])

		quake = {}
		quake['start'] = real_time[quake_start]
		quake['end'] = real_time[quake_end]

		local_slipping_times = velocity_norm[quake_start:quake_end, :] \
		> THRESHOLD_VELOCITY
		quake_displacement = displacement[quake_start:quake_end, :, :]
		slip = np.zeros(num_vertices)

		for i in range(num_vertices):
		if any(local_slipping_times[:, i]):
		starts = slip_beginnings(local_slipping_times[:, i])
		ends = slip_endings(local_slipping_times[:, i])
		slip[i] = np.linalg.norm(quake_displacement[ends, i, :]
		- quake_displacement[starts, i, :])

		quake['peakSlip'] = max(slip)
		quake['peakSlipRate'] = max(maximum_velocity[quake_start:quake_end])

		max_rupture_width = 0
		for time_step in range(local_slipping_times.shape[0]):
		slipping_time = local_slipping_times[time_step, :]
		if not any(slipping_time):
		continue

		slipping_coords = coordinates[slipping_time] \
		+ displacement[quake_start + time_step, slipping_time]
		if len(slipping_coords) > 1:
		pair = np.array(max_distance(slipping_coords))
		max_rupture_width = max(max_rupture_width,
		np.linalg.norm(pair[0] - pair[1]))

		quake['ruptureWidth'] = max_rupture_width
		quakes[quake_ID] = quake

		# output quake data to csv file
		csv_columns = quakes[0].keys()
		csv_file_path = os.path.join(out_path, 'events' + str(body_ID) + '.csv')
		try:
		with open(csv_file_path, 'w') as csv_file:
		writer = csv.DictWriter(csv_file, fieldnames=csv_columns)
		writer.writeheader()
		for quake_ID in quakes:
		writer.writerow(quakes[quake_ID])
		except IOError:
		print('I/O error')

		h5file.close()

data/tools/3d-performance.R

0 → 100644

+62 −0

Original line number	Diff line number	Diff line
		source('tools/support/findQuakes.R')

		finalTime <- 1000 # s
		convergenceVelocity <- 5e-5 # m/s

		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (basedir in c("rtol=1e-5_diam=1e-2")) {
		dir <- file.path(directories[['simulation']],
		'3d-lab', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
		relativeTime <- h5file['relativeTime'][]
		realTime <- finalTime * relativeTime

		velocityProxy<- h5file['/frictionalBoundary/velocity']

		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a very low velocity here.
		quakes <- findQuakes(1e-6 + convergenceVelocity, velocityProxy,
		## Note: We only look at the last 2000 timesteps here because
		## we're only interested in the last few events. This
		## dramatically reduces RAM usage and runtime.
		indices = seq(dim(velocityProxy)[1]-2000,
		dim(velocityProxy)[1]), c(1,3))
		quakes$beginning <- realTime[quakes$beginningIndex]
		quakes$ending <- realTime[quakes$endingIndex]
		quakes$duration <- quakes$ending - quakes$beginning
		numQuakes <- nrow(quakes)

		relaxedTime <- extendrange(c(quakes[[numQuakes-2,'beginning']],
		quakes[[numQuakes, 'ending']]), f=0.02)
		threeQuakeTimeMask <- (realTime > relaxedTime[[1]]) & (realTime < relaxedTime[[2]])
		plotMask <- threeQuakeTimeMask
		plotIndices<- which(plotMask)
		plotIndices<- c(plotIndices[1]-1,plotIndices,plotIndices[length(plotIndices)]+1)

		write(relaxedTime[[1]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'min', 'threequakes',
		basedir), 'tex')))
		write(relaxedTime[[2]],
		file.path(directories[['output']],
		paste.(pasteColon('timeframe', 'max', 'threequakes',
		basedir), 'tex')))

		timeWindow = realTime[plotIndices]

		relativeTau <- h5file['relativeTimeIncrement'][]
		fixedPointIterations <- h5file['/iterations/fixedPoint/final'][]
		multiGridIterations <- h5file['/iterations/multiGrid/final'][]
		write.csv(data.frame(time = timeWindow,
		timeIncrement = finalTime * relativeTau[plotIndices],
		fixedPointIterations = fixedPointIterations[plotIndices],
		multiGridIterations = multiGridIterations[plotIndices]),
		file.path(directories[['output']],
		paste.(pasteColon('3d-performance', basedir), 'csv')),
		row.names = FALSE, quote = FALSE)
		h5::h5close(h5file)
		}

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

0 → 100644

+67 −0

Original line number	Diff line number	Diff line
		source('tools/support/findQuakes.R')
		source('tools/support/writeContours.R')

		finalTime <- 1000 # s
		convergenceVelocity <- 5e-5 # m/s
		printlevels <- c('1','2','3','5',
		'10','20','30','50',
		'100','200','300','500',
		'1000')
		criticalVelocities <- 1e-6*as.numeric(printlevels) + convergenceVelocity

		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		directories <- ini::read.ini('config.ini')$directories
		dir.create(directories[['output']], recursive=TRUE, showWarnings=FALSE)
		for (basedir in c("rtol=1e-5_diam=1e-2")) {
		dir <- file.path(directories[['simulation']],
		'3d-lab', basedir)
		h5file <- h5::h5file(file.path(dir, 'output.h5'), 'r')
		relativeTime <- h5file['relativeTime'][]
		realTime <- finalTime * relativeTime

		velocityProxy<- h5file['/frictionalBoundary/velocity']

		## We are interested in an enlarged time range around actual events here,
		## (and no other quantities!) hence we pass a rather low velocity here.
		quakes <- findQuakes(200e-6 + convergenceVelocity, velocityProxy,
		## Note: We only look at the last 1000 timesteps here because
		## we're only interested in the last few events. This
		## dramatically reduces RAM usage and runtime.
		indices = seq(dim(velocityProxy)[1]-1000,
		dim(velocityProxy)[1]), c(1,3))
		quake <- quakes[nrow(quakes)-1,] # Q: why did we not need the -1 in julia?

		weakPatchGridVelocityProxy <- h5file["/weakPatchGrid/velocity"]

		stepSize <- 30 ## note: should/might differ by time/spatial resolution
		ts <- seq(quake$beginningIndex, quake$endingIndex, by=stepSize)
		dd <- data.frame(timeSteps = ts,
		times = realTime[ts],
		timeOffsets = realTime[ts] - realTime[quake$beginningIndex])

		fname = pasteColon('3d-velocity-contours', basedir)
		write.csv(dd, row.names = FALSE, quote = FALSE,
		file = file.path(directories[['output']],
		paste.(pasteColon(fname, 'times'), 'csv')))

		weakPatchGridXCoordinates <- h5file["/weakPatchGrid/xCoordinates"][]
		weakPatchGridZCoordinates <- h5file["/weakPatchGrid/zCoordinates"][]

		for (t in ts) {
		m <- weakPatchGridVelocityProxy[t,,,]
		s <- sqrt(m[,,,1]^2 + m[,,,3]^2)
		ret <- contourLines(weakPatchGridXCoordinates, weakPatchGridZCoordinates, s,
		level=criticalVelocities)

		for (i in seq(printlevels))
		writeContours(ret, criticalVelocities[[i]],
		file.path(directories[['output']],
		paste.(pasteColon(fname,
		'step', t,
		'level', printlevels[[i]]),
		'tex')))
		}
		h5::h5close(h5file)
		}

data/tools/comparison:lab-sim.R

0 → 100644

+42 −0

Original line number	Diff line number	Diff line
		paste. <- function(...) paste(..., sep='.')
		pasteColon<- function(...) paste(..., sep=':')

		basedir <- 'rfpitol=100e-7'

		directories <- ini::read.ini('config.ini')$directories
		labdata <- within(
		read.table(unz(file.path(directories[['experiment']],
		'B_Scale-model-earthquake-data.zip'),
		'scaleEQs_12-31_catalogue.txt'),
		sep='\t', header=FALSE, comment.char='%',
		col.names=c('frame',
		'meanSlipNature', 'meanSlipLab',
		'peakSlipNature', 'peakSlipLab',
		'ruptureWidthNature', 'ruptureWidthLab',
		'ignored', 'ignored', 'ignored', 'ignored', 'ignored')),
		{
		recurrenceLab <- c(NA, diff(frame))/10 # (10Hz camera: 10 frames ~ 1s)
		meanSlipLab <- meanSlipLab / 1e6 # micro m -> m
		peakSlipLab <- peakSlipLab / 1e6}) # micro m -> m
		## NB: We skip the first two quakes here, for compatibility with my old code
		## Maybe skipping the first was intentional so that each quake would have
		## a recurrence time. But skipping the second was certainly an accident
		labdata <- labdata[3:nrow(labdata),]

		simdata <- read.csv(file.path(directories[['output']],
		paste.(pasteColon('events', basedir), 'csv')))

		report <- function(lab, sim, quantity) {
		write.table(lab, file.path(directories[['output']],
		paste.(pasteColon('boxplot-data', 'lab',
		quantity), 'tex')),
		row.names=FALSE, col.names=FALSE)
		write.table(sim, file.path(directories[['output']],
		paste.(pasteColon('boxplot-data', 'simulation',
		basedir, quantity), 'tex')),
		row.names=FALSE, col.names=FALSE)
		}
		report(labdata$recurrenceLab, diff(simdata$beginning), 'recurrence')
		report(labdata$ruptureWidthLab, simdata$ruptureWidth, 'ruptureWidth')
		# meters -> millimeters
		report(labdata$peakSlipLab * 1000, simdata$peakSlip*1000, 'peakSlip')

data/tools/config.ini

0 → 100644

+9 −0

Original line number	Diff line number	Diff line
		#/home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/test/
		#/home/joscha/software/dune/build-debug/dune-tectonic/src/foam/output/test/
		#/home/joscha/Desktop/strikeslip/viscosity-1e3/
		#/home/joscha/Downloads/

		[directories]
		simulation = /home/joscha/software/dune/build-release/dune-tectonic/src/foam/output/tresca-0.6/
		experiment = ~/group/publications/2016-RosenauCorbiDominguezRudolfRitterPipping
		output = generated

data/tools/debug.py

0 → 100644

+135 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import numpy as np
		import csv
		import h5py
		import matplotlib.pyplot as plt

		from debug.outliers import outliers
		from debug.friction import truncated_friction
		from debug.state import aging_law
		from debug.diffplot import diffplot

		from support.maximum import maximum
		from support.norm import norm
		from support.find_quakes import find_quakes
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings
		from support.max_distance import max_distance

		from support.iterations import iterations

		NBODIES = 2
		FINAL_TIME = 15 # s
		FINAL_VELOCITY = 1e-5 # m/s
		THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY

		TANGENTIAL_COORDS = 1

		# friction params
		params = {
		'L' : 1e-5,
		'V0' : 1e-6,
		'mu0': 0.6,
		'a' : 0.010,
		'b' : 0.015
		}

		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('tools/config.ini')
		config.read(config_path)
		sim_path = config.get('directories', 'simulation')
		exp_path = config.get('directories', 'experiment')
		out_path = config.get('directories', 'output')

		# read hdf5 output file
		h5path = os.path.join(sim_path)
		h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		print(list(h5file.keys()))
		print(list(h5file['body1'].keys()))

		# read time
		relative_time = np.array(h5file['relativeTime'])
		relative_tau = np.array(h5file['relativeTimeIncrement'])
		relative_tau = np.delete(relative_tau, 0)

		real_time = relative_time * FINAL_TIME
		real_tau = relative_tau * FINAL_TIME

		print(len(relative_time))

		for body_ID in range(NBODIES):
		body = 'body' + str(body_ID)

		if body not in h5file:
		continue

		# velocity data
		v = abs(np.array(h5file[body + '/velocity']))

		# statistics
		avg_v = np.average(v[:,:,TANGENTIAL_COORDS], axis=1)
		min_v = np.min(v[:,:,TANGENTIAL_COORDS], axis=1)
		max_v = np.max(v[:,:,TANGENTIAL_COORDS], axis=1)

		# plot
		plt.figure(1)
		plt.subplot(311)
		plt.plot(min_v, color='gray', linestyle='--')
		plt.plot(avg_v, color='black', linestyle='-')
		plt.plot(max_v, color='gray', linestyle='--')
		plt.ylabel('slip rate')
		#-------------------------

		# state
		states = np.array(h5file[body + '/state'])
		states_calc = aging_law(params, states[0], v[:,:,TANGENTIAL_COORDS], real_tau)

		# statistics
		avg_states = np.average(states, axis=1)
		avg_states_calc = np.average(states_calc, axis=1)
		min_states = np.min(states, axis=1)
		max_states = np.max(states, axis=1)

		# plot
		plt.subplot(312)
		plt.plot(min_states, color='gray', linestyle='--')
		plt.plot(avg_states, color='black', linestyle='-')
		plt.plot(max_states, color='gray', linestyle='--')
		plt.plot(avg_states_calc, color='red', linestyle='-')
		plt.ylabel('state')
		#-------------------------

		# friction coefficient
		friction_coeff = np.array(h5file[body + '/coefficient'])

		#weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
		friction_coeff_calc = truncated_friction(params, v[:,:,TANGENTIAL_COORDS], states_calc)

		# statistics
		avg_friction_coeff = np.average(friction_coeff, axis=1)
		avg_friction_coeff_calc = np.average(friction_coeff_calc, axis=1)
		min_friction_coeff = np.min(friction_coeff, axis=1)
		max_friction_coeff = np.max(friction_coeff, axis=1)

		outliers_friction_coeff = outliers(avg_friction_coeff)

		# plot
		plt.subplot(313)
		plt.plot(min_friction_coeff, color='gray', linestyle='--')
		plt.plot(avg_friction_coeff, color='black', linestyle='-')
		plt.plot(max_friction_coeff, color='gray', linestyle='--')
		plt.plot(avg_friction_coeff_calc, color='red', linestyle='-')
		plt.plot(outliers_friction_coeff[0], outliers_friction_coeff[1], color='red', marker='+')
		plt.ylabel('friction coefficient')
		#-------------------------

		diffplot(friction_coeff[1], friction_coeff_calc[1], 'diff friction coeff')

		iterations(h5file, FINAL_TIME)

		plt.show()

		h5file.close()
		No newline at end of file

data/tools/debug/init.py

0 → 100644

+0 −0

Original line number	Diff line number	Diff line

data/tools/debug/diffplot.py

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		import numpy as np
		import matplotlib.pyplot as plt

		def diffplot(v1, v2, title):
		v_diff = v1 - v2
		t = np.linspace(0, len(v_diff)-1, len(v_diff), endpoint=True)

		# plot
		fig = plt.figure()
		ax = fig.add_subplot(1, 1, 1)
		ax.plot(t, v1, color='black', linestyle='--')
		ax.plot(t, v2, color='gray', linestyle='--')
		ax.plot(t, v_diff, color='red', linestyle='-')
		ax.set_ylabel(title)
		#-------------------------

		fig.canvas.draw()
		No newline at end of file

data/tools/debug/friction.py

0 → 100644

+8 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def truncated_friction(params, v, alpha):
		vmin = params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )
		clipped_v = (v/vmin).clip(1)
		return params['a'] * np.log(clipped_v)

		#return (params['a'] * np.log(v/params['V0']) - params['mu0'] - (params['b']/params['a'])*np.array(alpha) ).clip(0)
		No newline at end of file

data/tools/debug/outliers.py

0 → 100644

+6 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def outliers(data):
		index = [i for i,x in enumerate(data) if np.abs(x)==np.inf]
		val = np.zeros(len(index))
		return [index, val]
		No newline at end of file

data/tools/debug/state.py

0 → 100644

+25 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def lift_singularity(c, x):
		def local_lift(y):
		if (y <= 0):
		return -c
		else:
		return -np.expm1(c * y) / y

		return np.array([local_lift(y) for y in x])


		def aging_law(params, initial_alpha, v, time_steps):
		#auto tangentVelocity = velocity_field[localToGlobal_[i]];
		#tangentVelocity[0] = 0.0;
		#double const V = tangentVelocity.two_norm();

		alpha = [initial_alpha]

		for i,tau in enumerate(time_steps):
		mtoL = -tau / params['L']
		next_alpha = np.log(np.exp(alpha[-1] + v[i] * mtoL) + params['V0'] * lift_singularity(mtoL, v[i]))
		alpha.append(next_alpha)

		return alpha
		No newline at end of file

data/tools/elias.py

0 → 100644

+73 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import numpy as np
		import csv
		import h5py
		import matplotlib.pyplot as plt

		from debug.outliers import outliers
		from debug.friction import truncated_friction
		from debug.state import aging_law
		from debug.diffplot import diffplot

		from support.maximum import maximum
		from support.norm import norm
		from support.find_quakes import find_quakes
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings
		from support.max_distance import max_distance

		from support.iterations import iterations
		from support.friction_stats import friction_stats

		def build_patch(coords, percentage):
		x_coords = coords[:, 0]
		xmin = np.min(x_coords)
		xmax = np.max(x_coords)
		delta_x = (1 - percentage)*(xmax - xmin)/2

		xmin = xmin + delta_x
		xmax = xmax - delta_x

		return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]

		FINAL_TIME = 1000 # s
		FINAL_VELOCITY = 5e-5 # m/s
		THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY

		TANGENTIAL_COORDS = 0

		# friction params
		params = {
		'L' : 1e-5,
		'V0' : 1e-6,
		'mu0': 0.6,
		'a' : 0.010,
		'b' : 0.015
		}

		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('tools/config.ini')
		config.read(config_path)
		sim_path = config.get('directories', 'simulation')
		exp_path = config.get('directories', 'experiment')
		out_path = config.get('directories', 'output')

		# read hdf5 output file
		h5path = os.path.join(sim_path)
		h5file = h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		print(list(h5file.keys()))
		print(list(h5file['frictionalBoundary'].keys()))

		iterations(h5file, FINAL_TIME)

		coords = np.array(h5file['frictionalBoundary/coordinates'])
		patch = build_patch(coords, 0.05)

		friction_stats(h5file, 0, FINAL_TIME, [], [0, 50], TANGENTIAL_COORDS)

		plt.show()

		h5file.close()
		No newline at end of file

data/tools/generate-2d.bash

0 → 100644

+25 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env bash

		set -e

		rr() {
		echo "$(date)" Running: Rscript $@
		Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
		}

		# contours
		rr 2d-velocity-contours.R

		# dip
		rr 2d-dip-contours.R

		# iterations / adaptivity
		rr 2d-performance.R

		# box plot (one half)
		rr 2d-event-writer.R

		# fpi data
		rr 2d-fpi-tolerance.R

		date

data/tools/generate-3d.bash

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env bash

		set -e

		rr() {
		echo "$(date)" Running: Rscript $@
		Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
		}

		# performance
		rr 3d-performance.R

		# velocity contours
		rr 3d-velocity-contours.R

		date

data/tools/generate-others.bash

0 → 100644

+13 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env bash

		set -e

		rr() {
		echo "$(date)" Running: Rscript $@
		Rscript --vanilla --default-packages=grDevices,methods,stats,utils tools/$@
		}

		# box plot (one half)
		rr comparison:lab-sim.R

		date

data/tools/generated/events1.csv

0 → 100644

+20 −0

Original line number	Diff line number	Diff line
		start,peakSlipRate,end,ruptureWidth,peakSlip
		20.086059570312422,5.7045853762396602e-06,20.086669921874922,0,0.0
		20.090332031249918,5.6219949973608878e-06,20.092163085937418,0,0.0
		20.094604492187415,5.8412336903395239e-06,20.095825195312415,0,0.0
		20.099487304687411,5.3352579051826403e-06,20.100097656249908,0,0.0
		20.143737792968619,5.4285740142335498e-06,20.144958496093619,0,0.0
		20.148620605468615,5.6370561043864526e-06,20.149536132812365,0,0.0
		20.353851318359059,6.0525078256123743e-06,20.354461669921559,0,0.0
		20.358734130859055,5.1806538408711459e-06,20.359344482421555,0,0.0
		20.371856689452795,5.6547173619088325e-06,20.373077392577795,0,0.0
		20.425262451171495,6.4545930360194781e-06,20.428314208983991,0.088385894124722503,1.4298690605481336e-08
		20.479583740233945,5.8310594282664788e-06,20.481414794921445,0,0.0
		20.483245849608942,5.2537138171449348e-06,20.485076904296442,0,6.2000101383992609e-09
		20.531158447265149,5.2316818926703734e-06,20.531768798827649,0,0.0
		20.544586181640138,5.4092184353134343e-06,20.545196533202635,0,0.0
		20.593109130858842,6.190190600109647e-06,20.595550537108842,0,0.0
		20.608673095702581,5.5618674305029011e-06,20.609283447265081,0,0.0
		20.625152587890067,5.4920875934131527e-06,20.625762939452567,0,0.0
		20.630035400390064,5.1740512343599518e-06,20.630645751952564,0,0.0
		20.63369750976506,5.7392806837644885e-06,20.63491821289006,0,0.0

data/tools/main.py

0 → 100644

+71 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import numpy as np
		import csv
		import h5py
		import matplotlib.pyplot as plt

		from debug.outliers import outliers
		from debug.friction import truncated_friction
		from debug.state import aging_law
		from debug.diffplot import diffplot

		from support.maximum import maximum
		from support.norm import norm
		from support.find_quakes import find_quakes
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings
		from support.max_distance import max_distance

		from support.io import read_h5file

		from support.iterations import iterations
		from support.friction_stats import friction_stats

		def build_patch(coords, percentage):
		x_coords = coords[:, 0]
		xmin = np.min(x_coords)
		xmax = np.max(x_coords)
		delta_x = (1 - percentage)*(xmax - xmin)/2

		xmin = xmin + delta_x
		xmax = xmax - delta_x

		return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]


		NBODIES = 2
		FINAL_TIME = 150 # s
		FINAL_VELOCITY = 2e-4 # m/s
		THRESHOLD_VELOCITY = 0.5*FINAL_VELOCITY # 1000e-6 + FINAL_VELOCITY

		TANGENTIAL_COORDS = 1

		# friction params
		params = {
		'L' : 1e-5,
		'V0' : 1e-6,
		'mu0': 0.6,
		'a' : 0.010,
		'b' : 0.015
		}

		h5file = read_h5file()
		print(list(h5file.keys()))

		iterations(h5file, FINAL_TIME)

		for body_ID in range(NBODIES):
		body = 'body' + str(body_ID)

		if body not in h5file:
		continue

		coords = np.array(h5file[body + '/coordinates'])
		patch = build_patch(coords, 0.05)

		friction_stats(h5file, body_ID, FINAL_TIME, [], [0, FINAL_TIME])

		plt.show()

		h5file.close()
		No newline at end of file

data/tools/support/init.py

0 → 100644

+1 −0

Original line number	Diff line number	Diff line

data/tools/support/init.pyc

0 → 100644

+106 B

File added.

No diff preview for this file type.

View file

data/tools/support/find_quakes.py

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings


		def find_quakes(threshold_velocity, maximum_velocities):
		slipping_times = maximum_velocities > threshold_velocity

		quake_starts = slip_beginnings(slipping_times)
		quake_ends = slip_endings(slipping_times)

		# remove incomplete quakes
		min_len = min(len(quake_starts), len(quake_ends))
		quake_ends = quake_ends[0:min_len]
		quake_starts = quake_starts[0:min_len]

		return [quake_starts, quake_ends]

data/tools/support/find_quakes.pyc

0 → 100644

+741 B

File added.

No diff preview for this file type.

View file

data/tools/support/find_quakes.py~

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings


		def find_quakes(threshold_velocity, maximum_velocities):
		slipping_times = maximum_velocities > threshold_velocity
		print(slipping_times)

		quake_starts = slip_beginnings(slipping_times)
		quake_ends = slip_endings(slipping_times)

		# remove incomplete quakes
		min_len = min(len(quake_starts), len(quake_ends))
		quake_ends = quake_ends[0:min_len]
		quake_starts = quake_starts[0:min_len]

		return [quake_starts, quake_ends]

data/tools/support/friction_stats.py

0 → 100644

+74 −0

Original line number	Diff line number	Diff line
		import numpy as np
		import matplotlib.pyplot as plt

		def friction_stats(h5file, body_ID, FINAL_TIME, patch=[], interval=[], TANGENTIAL_COORDS=1):
		body = 'body' + str(body_ID) # 'frictionalBoundary' 'body' + str(body_ID)

		coords = np.array(h5file[body + '/coordinates'])
		if len(patch) == 0:
		patch = np.linspace(0, len(coords)-1, len(coords), endpoint=True, dtype='int8')

		# read time
		time = np.array(h5file['relativeTime']) * FINAL_TIME
		time = np.delete(time, 0)
		if len(interval) == 0:
		interval = [0, FINAL_TIME]
		t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
		time = time[t]

		fig = plt.figure()

		# velocity data
		v = abs(np.array(h5file[body + '/velocity']))
		v_t = v[t,:,TANGENTIAL_COORDS]
		v_tx = v_t[:,patch]
		# statistics
		avg_v = np.average(v_tx, axis=1)
		min_v = np.min(v_tx, axis=1)
		max_v = np.max(v_tx, axis=1)
		# plot
		ax_slip = fig.add_subplot(3, 1, 1)
		#ax_slip.plot(time, min_v, color='gray', linestyle='--')
		ax_slip.plot(time, avg_v, color='black', linestyle='-')
		#ax_slip.plot(time, max_v, color='gray', linestyle='--')
		ax_slip.set_ylabel('slip rate')
		#ax_slip.set_yscale('log')
		#-------------------------

		print(np.min(min_v))
		print(np.average(avg_v))
		print(np.max(max_v))

		# state
		states = np.array(h5file[body + '/state'])
		states_t = states[t,:]
		states_tx = states_t[:,patch]
		# statistics
		avg_states = np.average(states_tx, axis=1)
		min_states = np.min(states_tx, axis=1)
		max_states = np.max(states_tx, axis=1)
		# plot
		ax_state = fig.add_subplot(3, 1, 2)
		ax_state.plot(time, min_states, color='gray', linestyle='--')
		ax_state.plot(time, avg_states, color='black', linestyle='-')
		ax_state.plot(time, max_states, color='gray', linestyle='--')
		ax_state.set_ylabel('state')
		#-------------------------

		# friction coefficient
		friction_coeff = np.array(h5file[body + '/coefficient'])
		friction_coeff_t = friction_coeff[t,:]
		friction_coeff_tx = friction_coeff_t[:,patch]
		# statistics
		avg_friction_coeff = np.average(friction_coeff_tx, axis=1)
		min_friction_coeff = np.min(friction_coeff_tx, axis=1)
		max_friction_coeff = np.max(friction_coeff_tx, axis=1)
		# plot
		ax_friction = fig.add_subplot(3, 1, 3)
		ax_friction.plot(time, min_friction_coeff, color='gray', linestyle='--')
		ax_friction.plot(time, avg_friction_coeff, color='black', linestyle='-')
		ax_friction.plot(time, max_friction_coeff, color='gray', linestyle='--')
		ax_friction.set_ylabel('friction coefficient')
		#-------------------------

		fig.canvas.draw()
		No newline at end of file

data/tools/support/io.py

0 → 100644

+49 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import h5py

		def remove_comment(str, char = "#"):
		split_str = str.split(char, 1)
		return split_str[0]

		# tag = 'simulation', 'experiment', 'output'
		def read_h5file(tag = 'simulation'):
		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('tools/config.ini')
		config.read(config_path)
		path = config.get('directories', tag)

		# read hdf5 output file
		h5path = os.path.join(path)

		return h5py.File(os.path.join(h5path, 'output.h5'), 'r')

		# tag = 'simulation', 'experiment', 'output'
		def read_params(file_name, tag = 'simulation'):
		# read config ini
		config = cp.ConfigParser()
		config_path = os.path.join('tools/config.ini')
		config.read(config_path)
		path = config.get('directories', tag)

		# read cfg parameter file
		params = cp.ConfigParser()
		params_path = os.path.join(path, file_name)
		params.read(params_path)

		#with open(params_path) as stream:
		# params.read_string("[top]\n" + stream.read())

		out = {
		'L' : float(remove_comment(params.get('boundary.friction', 'L'))),
		'V0' : float(remove_comment(params.get('boundary.friction', 'V0'))),
		'mu0': float(remove_comment(params.get('boundary.friction', 'mu0'))),
		'a' : float(remove_comment(params.get('boundary.friction.weakening', 'a'))),
		'b' : float(remove_comment(params.get('boundary.friction.weakening', 'b'))),
		'bodyCount' : int(remove_comment(params.get('problem', 'bodyCount'))),
		'finalTime' : float(remove_comment(params.get('problem', 'finalTime'))),
		'finalVelocity' : float(remove_comment(params.get('boundary.dirichlet', 'finalVelocity')))
		}

		return out
		No newline at end of file

data/tools/support/iterations.py

0 → 100644

+45 −0

Original line number	Diff line number	Diff line
		import numpy as np
		import matplotlib.pyplot as plt

		def iterations(h5file, FINAL_TIME, interval = []):

		# read time
		time = np.array(h5file['relativeTime']) * FINAL_TIME
		time = np.delete(time, 0)
		if len(interval) == 0:
		interval = [0, FINAL_TIME]
		t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
		time = time[t]

		tau = np.array(h5file['relativeTimeIncrement']) * FINAL_TIME
		tau = np.delete(tau, 0)

		# read fpi iterations
		fpi_final = np.array(h5file['iterations/fixedPoint/final'])
		fpi_final = np.delete(fpi_final, 0)
		#fpi_total = np.array(h5file['iterations/fixedPoint/total'])

		# read multigrid iterations
		multigrid_final = np.array(h5file['iterations/multiGrid/final'])
		multigrid_final = np.delete(multigrid_final, 0)
		#multigrid_total = np.array(h5file['iterations/multiGrid/total'])

		# plot
		fig = plt.figure()

		ax_fpi = fig.add_subplot(3, 1, 1)
		ax_fpi.plot(time, fpi_final[t], color='black', linestyle='-')
		ax_fpi.set_ylabel('fpi')
		#-------------------------

		ax_mg = fig.add_subplot(3, 1, 2)
		ax_mg.plot(time, multigrid_final[t], color='black', linestyle='-')
		ax_mg.set_ylabel('multigrid')
		#-------------------------

		ax_tau = fig.add_subplot(3, 1, 3)
		ax_tau.plot(time, tau[t], color='black', linestyle='-')
		ax_tau.set_ylabel('tau')
		#-------------------------

		fig.canvas.draw()
		No newline at end of file

data/tools/support/max_distance.py

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		from itertools import combinations


		def square_distance(x, y):
		return sum([(xi-yi)**2 for xi, yi in zip(x, y)])


		def max_distance(points):
		max_square_distance = 0
		max_pair = (0, 0)
		for pair in combinations(points, 2):
		sq_dist = square_distance(*pair)

		if sq_dist > max_square_distance:
		max_square_distance = sq_dist
		max_pair = pair
		return max_pair

data/tools/support/max_distance.pyc

0 → 100644

+884 B

File added.

No diff preview for this file type.

View file

data/tools/support/max_distance.py~

0 → 100644

+17 −0

Original line number	Diff line number	Diff line
		from itertools import combinations


		def square_distance(x, y):
		return sum([(xi-yi)**2 for xi, yi in zip(x, y)])


		def max_distance(points):
		max_square_distance = 0
		max_pair = (0, 0)
		for pair in combinations(points, 2):
		sq_dist = square_distance(*pair)
		print(sq_dist)
		if sq_dist > max_square_distance:
		max_square_distance = sq_dist
		max_pair = pair
		return max_pair

data/tools/support/max_velocity.py

0 → 100644

+4 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def max_velocity(x) :
		return maximum(norm(x))
		No newline at end of file

data/tools/support/max_velocity.pyc

0 → 100644

+649 B

File added.

No diff preview for this file type.

View file

data/tools/support/max_velocity.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def max_velocity(x) :
		size = x.shape
		res = np.zeros(size[0])

		for time_step in range(size[0]) :
		for vertex in range(size[1]) :
		res[time_step] = max(res[time_step], np.linalg.norm(x[time_step, vertex, :]))

		return res
		No newline at end of file

data/tools/support/maximum.py

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np


		def maximum(x):
		size = x.shape
		res = np.zeros(size[0])

		for time_step in range(size[0]):
		res[time_step] = max(res[time_step], max(x[time_step, :]))

		return res

data/tools/support/maximum.pyc

0 → 100644

+560 B

File added.

No diff preview for this file type.

View file

data/tools/support/maximum.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def maximum(x) :
		size = x.shape
		res = np.zeros(size[0])

		for time_step in range(size[0]) :
		for vertex in range(size[1]) :
		res[time_step] = max(res[time_step], max(x[time_step, vertex]))

		return res
		No newline at end of file

data/tools/support/norm.py

0 → 100644

+12 −0

Original line number	Diff line number	Diff line
		import numpy as np


		def norm(x):
		size = x.shape
		res = np.zeros((size[0], size[1]))

		for time_step in range(size[0]):
		for vertex in range(size[1]):
		res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])

		return res

data/tools/support/norm.pyc

0 → 100644

+616 B

File added.

No diff preview for this file type.

View file

data/tools/support/norm.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

		def norm(x) :
		size = x.shape
		res = np.zeros(size[0], size[1])

		for time_step in range(size[0]) :
		for vertex in range(size[1]) :
		res[time_step, vertex] = np.linalg.norm(x[time_step, vertex, :])

		return res

data/tools/support/slip_beginnings.py

0 → 100644

+15 −0

Original line number	Diff line number	Diff line
		import array as ar


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

		length = 0
		for i in range(len(x)):
		if (not prev and x[i]):
		starts[length] = i
		length += 1
		prev = x[i]
		return starts[:length]

data/tools/support/slip_beginnings.pyc

0 → 100644

+837 B

File added.

No diff preview for this file type.

View file

data/tools/support/slip_beginnings.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

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

data/tools/support/slip_endings.py

0 → 100644

+16 −0

Original line number	Diff line number	Diff line
		import array as ar


		def slip_endings(x):
		# returns indicies i for which x[i-1]>0 and x[i]=0
		ends = ar.array('i', (0 for i in range(len(x))))
		prev = False

		length = 0
		for i in range(len(x)):
		if (prev and not x[i]):
		ends[length] = i
		length += 1
		prev = x[i]

		return ends[:length]

data/tools/support/slip_endings.pyc

0 → 100644

+823 B

File added.

No diff preview for this file type.

View file

data/tools/support/slip_endings.py~

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		import numpy as np

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

data/tools/support/trapezoidal.py

0 → 100644

+8 −0

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

data/tools/support/trapezoidal.py~

0 → 100644

+8 −0

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

data/tools/support/write_contours.py

0 → 100644

+11 −0

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

data/tools/truncation.py

0 → 100644

+116 −0

Original line number	Diff line number	Diff line
		import configparser as cp
		import os
		import numpy as np
		import csv
		import h5py
		import matplotlib.pyplot as plt

		from debug.outliers import outliers
		from debug.friction import truncated_friction
		from debug.state import aging_law
		from debug.diffplot import diffplot

		from support.maximum import maximum
		from support.norm import norm
		from support.find_quakes import find_quakes
		from support.slip_beginnings import slip_beginnings
		from support.slip_endings import slip_endings
		from support.max_distance import max_distance

		from support.io import read_h5file
		from support.io import read_params

		from support.iterations import iterations
		from support.friction_stats import friction_stats

		def build_time(h5file, final_time, interval = []):
		# read time
		time = np.array(h5file['relativeTime']) * FINAL_TIME
		time = np.delete(time, 0)

		if len(interval) == 0:
		interval = [0, FINAL_TIME]

		t = [i for i in range(len(time)) if time[i]>=interval[0] and time[i]<=interval[1]]
		return t, time[t]

		def build_patch(coords, percentage = 1.0):
		x_coords = coords[:, 0]
		xmin = np.min(x_coords)
		xmax = np.max(x_coords)
		delta_x = (1 - percentage)*(xmax - xmin)/2

		xmin = xmin + delta_x
		xmax = xmax - delta_x

		return [i for i in range(len(x_coords)) if x_coords[i]>=xmin and x_coords[i]<=xmax]

		# read h5 output file
		h5file = read_h5file()
		print(list(h5file.keys()))

		# read problem parameters
		params = read_params('foam.cfg')
		print(params)

		TANGENTIAL_COORDS = 1
		FINAL_TIME = params['finalTime']
		NBODIES = params['bodyCount']

		t, time = build_time(h5file, FINAL_TIME, [0, FINAL_TIME])

		for body_ID in range(NBODIES):
		body = 'body' + str(body_ID)

		if body not in h5file:
		continue

		coords = np.array(h5file[body + '/coordinates'])
		patch = build_patch(coords, 1.0)

		fig = plt.figure()

		# velocity and state data
		v = abs(np.array(h5file[body + '/velocity']))[:,:,TANGENTIAL_COORDS]
		alpha = np.array(h5file[body + '/state'])

		vmin = 0 #params['V0'] / np.exp( ( params['mu0'] + params['b'] * np.array(alpha)) / params['a'] )

		v_truncation = (v<=vmin)[t,:]
		v_truncated = np.sum(v_truncation, axis=1)

		v_comp = (np.abs(v-vmin)<1e-14)[t,:]
		v_small = np.sum(v_comp, axis=1)

		weighted_normal_stress = np.array(h5file[body + '/weightedNormalStress'])
		second_deriv_large = (v<0)[t,:] #((- weighted_normal_stress[None, ...] * params['a'] / v)>1e10)[0,t,:]
		second_deriv = np.sum(second_deriv_large, axis=1)

		total_truncated = np.sum(v_truncation \| v_comp \| second_deriv_large, axis=1)

		# plot v_small
		ax_v_small = fig.add_subplot(4, 1, 1)
		ax_v_small.plot(time, v_small, color='black', linestyle='-')
		ax_v_small.set_ylabel('# v-vmin < 1e-14')
		#-------------------------

		# plot v_truncated
		ax_v_truncated = fig.add_subplot(4, 1, 2 )
		ax_v_truncated.plot(time, v_truncated, color='black', linestyle='-')
		ax_v_truncated.set_ylabel('# v <= vmin')
		#-------------------------

		# plot second_deriv_large
		ax_second_deriv = fig.add_subplot(4, 1, 3)
		ax_second_deriv.plot(time, second_deriv, color='black', linestyle='-')
		ax_second_deriv.set_ylabel('# second_deriv large')
		#-------------------------

		# plot total_truncated
		ax_truncated = fig.add_subplot(4, 1, 4)
		ax_truncated.plot(time, total_truncated, color='black', linestyle='-')
		ax_truncated.set_ylabel('# total truncated')
		#-------------------------
		plt.show()

		h5file.close()
		No newline at end of file

debugging.m

0 → 100644

+589 −0

Original line number	Diff line number	Diff line
		% householder reflection

		v0 = [1; 0];
		v0 = v0./norm(v0);

		v1 = [1; -1];
		v1 = v1./norm(v1);

		v = 1/2.*(v1-v0);
		vnorm = norm(v);

		H = eye(2) - 2/(vnorm^2).(vv')

		A = [...
		];

		f = [...
		0;
		-4.33681e-19;
		0;
		1.0842e-19;
		-2.66671e-06;
		-4.94914e-07;
		0;
		-2.71051e-20;
		-4.19923e-06;
		-4.21374e-07;
		-3.92235e-06;
		-8.64817e-07;
		0;
		-5.42101e-20;
		-4.96301e-06;
		-4.13502e-06;
		0;
		0;
		-7.03515e-06;
		-5.84881e-06;
		-4.93309e-06;
		-2.37804e-06;
		0;
		0;
		-2.90755e-06;
		-1.23334e-06;
		0;
		0;
		-2.2449e-06;
		-8.21563e-07;
		0;
		0;
		0;
		0;
		0;
		0;
		-7.24652e-06;
		-8.25293e-07;
		0;
		2.4104e-06;
		-8.19294e-06;
		-1.18617e-06;
		-5.58223e-06;
		-1.69348e-06;
		0;
		1.88001e-06;
		-1.03387e-05;
		-8.22886e-06;
		0;
		0;
		-1.07258e-05;
		-8.62344e-06;
		-9.76077e-06;
		-4.93954e-06;
		-9.10261e-06;
		-7.35307e-06;
		-4.01973e-06;
		-2.38351e-06;
		-5.92659e-06;
		-4.23858e-07;
		-7.9911e-06;
		-4.04139e-06;
		-3.80594e-06;
		-3.06903e-06;
		-3.92924e-06;
		-2.71681e-06;
		-5.12975e-06;
		-1.37642e-06;
		-4.79892e-06;
		-1.14935e-06;
		-7.19726e-06;
		-3.54967e-06;
		-6.62755e-06;
		-6.1499e-07;
		0;
		1.28701e-05;
		0;
		0;
		-1.09944e-05;
		-1.92221e-06;
		0;
		1.17853e-05;
		-1.21456e-05;
		-2.20784e-06;
		-8.3682e-06;
		-2.01449e-06;
		0;
		0;
		-1.43122e-05;
		-1.00938e-05;
		0;
		0;
		-1.48374e-05;
		-1.03947e-05;
		-1.3684e-05;
		-6.31825e-06;
		-1.35318e-05;
		-9.3703e-06;
		-5.76679e-06;
		-1.88878e-06;
		-9.18524e-06;
		-1.51544e-06;
		-1.16296e-05;
		-5.64398e-06;
		-5.20301e-06;
		-2.13559e-06;
		-5.55432e-06;
		-2.017e-06;
		-7.96841e-06;
		-1.7895e-06;
		-7.95375e-06;
		-1.76175e-06;
		-1.0596e-05;
		-5.24553e-06;
		-1.02345e-05;
		-1.74422e-06;
		-2.16064e-05;
		-1.85059e-06;
		-2.24752e-05;
		-3.34957e-07;
		-1.34774e-05;
		-2.51353e-06;
		-2.18312e-05;
		-1.04202e-06;
		-1.99898e-05;
		-1.463e-06;
		-1.94256e-05;
		-2.11821e-06;
		-5e-05;
		5.03269e-07;
		-2.00631e-05;
		-2.41829e-06;
		-2.11637e-05;
		-1.12298e-06;
		-1.61508e-05;
		-1.11752e-05;
		-1.7267e-05;
		-6.9337e-06;
		-1.43931e-05;
		-6.46903e-06;
		-1.70354e-05;
		-2.45177e-06;
		-5e-05;
		-2.50465e-06;
		-5e-05;
		-1.21588e-05;
		-2.01431e-05;
		-1.19797e-05;
		];

		d = 0.00075745;
		A2 = [...
		];

		f2 = [...
		0;
		0;
		0;
		0;
		-3.01994e-06;
		-1.76298e-07;
		0;
		0;
		-4.78204e-06;
		-2.0995e-07;
		-5.01951e-06;
		1.22653e-06;
		0;
		0;
		-4.98476e-06;
		-3.57931e-06;
		0;
		0;
		-7.06572e-06;
		-4.97786e-06;
		-5.40025e-06;
		-1.91591e-06;
		0;
		0;
		-3.26344e-06;
		-9.4196e-07;
		0;
		0;
		-3.47107e-06;
		1.69693e-06;
		0;
		0;
		0;
		0;
		0;
		0;
		-8.36624e-06;
		-3.21786e-07;
		0;
		0;
		-9.19771e-06;
		-4.8646e-07;
		-9.46319e-06;
		2.13218e-06;
		0;
		0;
		-1.03261e-05;
		-6.78209e-06;
		0;
		0;
		-1.09634e-05;
		-7.09926e-06;
		-1.01615e-05;
		-3.77678e-06;
		-9.09151e-06;
		-6.13484e-06;
		-9.06507e-06;
		4.28268e-06;
		-6.83374e-06;
		-2.4706e-07;
		-8.51881e-06;
		-3.22169e-06;
		-7.83369e-06;
		3.49294e-06;
		-8.45665e-06;
		3.90404e-06;
		-7.913e-06;
		1.91502e-06;
		-7.13457e-06;
		1.77604e-06;
		-7.71821e-06;
		-2.87949e-06;
		-7.58455e-06;
		-3.19117e-07;
		0;
		0;
		0;
		0;
		-1.17905e-05;
		-6.57623e-07;
		0;
		0;
		-1.26738e-05;
		-7.75536e-07;
		-1.28663e-05;
		2.59111e-06;
		0;
		0;
		-1.33726e-05;
		-8.23012e-06;
		0;
		0;
		-1.40382e-05;
		-8.4827e-06;
		-1.3445e-05;
		-4.60248e-06;
		-1.21277e-05;
		-7.7075e-06;
		-1.20933e-05;
		5.59837e-06;
		-1.00511e-05;
		-5.83173e-07;
		-1.16989e-05;
		-4.19969e-06;
		-1.07938e-05;
		5.03904e-06;
		-1.14551e-05;
		5.33497e-06;
		-1.1057e-05;
		2.39173e-06;
		-1.01977e-05;
		2.22842e-06;
		-1.08385e-05;
		-3.93909e-06;
		-1.09119e-05;
		-6.59537e-07;
		-1.51572e-05;
		6.80143e-06;
		-1.97651e-05;
		7.92292e-06;
		-1.3582e-05;
		-8.81729e-07;
		-1.73372e-05;
		7.52177e-06;
		-1.63794e-05;
		2.79237e-06;
		-1.33659e-05;
		2.38778e-06;
		-5e-05;
		7.93536e-06;
		-1.92072e-05;
		-9.96979e-07;
		-1.93357e-05;
		3.15138e-06;
		-1.54595e-05;
		-9.11795e-06;
		-1.68142e-05;
		-5.1397e-06;
		-1.4113e-05;
		-4.66534e-06;
		-1.63265e-05;
		-9.49749e-07;
		-5e-05;
		-1.0631e-06;
		-5e-05;
		-1.02021e-05;
		-1.98778e-05;
		-9.99727e-06;
		];

		d2 = 0;

		sum(sum(abs(f-f2)>1e-14))


		sum(sum(abs(A-A2)>1e-14))


		alpha0 = 100;
		V = 1e-5; %1e-20;
		V0 = 5e-5;
		L = 2.25e-5;

		alpha = @(t) log(exp(alpha0 - V.t./L ) - V0 (exp(-t./L.*V)-1)./V);

		x = [0:0.001:0.1];
		alpha(x)
		plot(x, alpha(x));

		%abs(d-d2)

		% x = [...
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% -1.2337e-10;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% -1.2337e-10;
		% 0;
		% -1.2337e-10;
		% 0;
		% 0;
		% 0;
		% ];
		%
		% ignore = [...
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 1;
		% 1;
		% 0;
		% 0;
		% 1;
		% 1;
		% 0;
		% 0;
		% 1;
		% 1;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 1;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 0;
		% 1;
		% 0;
		% 1;
		% 0;
		% 0;
		% 0;
		% ];
		%
		% newA = A;
		% newf = f;
		% for i=1:length(ignore)
		% if (ignore(i)==1)
		% for j=1:length(ignore)
		% newA(i,j) = (i==j);
		% end
		% newf(i) = x(i);
		% end
		% end
		% newA
		% newf

		% ignore = zeros(1, length(f));
		% for i=1:length(ignore)
		% if (A(i,i)==0)
		% ignore(1, i) = 1;
		% end
		% end
		%
		% indices = [];
		% for i=1:length(ignore)
		% if (ignore(i)==0)
		% indices = [indices i];
		% end
		% end
		%
		% len = length(indices);
		% newA = zeros(len, len);
		% newf = zeros(1, len);
		% for i=1:len
		% newf(i) = f(indices(i));
		%
		% for j=1:len
		% newA(i,j) = A(indices(i), indices(j));
		% end
		% end
		% newf;
		% newA

		%inv(A);

		% sol = newA\newf
		%
		% tnnmgSol = [...
		% -1.03349;
		% -0.0502635;
		% 0.0161336;
		% -0.00852914;
		% -0.0130659;
		% 0.00660115;
		% -0.434915;
		% -0.0128241;
		% -0.00893317;
		% 0.00884256;
		% -0.0110232;
		% -0.10936;
		% 1.14978;
		% 0.0581132;
		% 0.019598;
		% 0.190758;
		% 0.478393;
		% -0.00433049;
		% 0.0513446;
		% 0.353218;
		% -0.00482644;
		% 0.11758;
		% 0;
		% 0;
		% -0.013677;
		% 0.0512563;
		% 0;
		% 0;
		% -0.00514197;
		% -0.117664;
		% 0;
		% 0;
		% 0.130034;
		% 0.574041;
		% 0.0136965;
		% 0.214953;
		% 0.0335986;
		% -0.00703477;
		% 0.0438202;
		% 0.346506;
		% -0.000896805;
		% 0.127837;
		% 0.041946;
		% 0.217409;
		% -1.2337e-10;
		% 0.112528;
		% -0.00828618;
		% 0.00453345;
		% -0.00898709;
		% 0.0781683;
		% 0.092355;
		% -0.556143;
		% -0.00443686;
		% -0.140953;
		% 0.0424371;
		% -0.250082;
		% -0.00336878;
		% 0.00105471;
		% -1.2337e-10;
		% 0.0188452;
		% -1.2337e-10;
		% -0.14111;
		% 0.00298499;
		% -0.190904;
		% ];
		%
		% sol-tnnmgSol
		No newline at end of file

dune-tectonic.pc.in

+1 −1

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

dune.module

+1 −1

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

dune/tectonic/CMakeLists.txt

+24 −14

Original line number	Diff line number	Diff line
		add_subdirectory("data-structures")
		add_subdirectory("factories")
		add_subdirectory("io")
		add_subdirectory("multi-threading")
		add_subdirectory("problem-data")
		add_subdirectory("spatial-solving")
		add_subdirectory("tests")
		add_subdirectory("time-stepping")
		add_subdirectory("utils")

		add_custom_target(tectonic_dune SOURCES
		assemblers.hh
		assemblers.cc
		explicitgrid.hh
		explicitvectors.hh
		gridselector.hh
		)

		#install headers
		install(FILES
		body.hh
		frictiondata.hh
		frictionpotential.hh
		globalfrictiondata.hh
		globalfriction.hh
		globalratestatefriction.hh
		gravity.hh
		localfriction.hh
		minimisation.hh
		myblockproblem.hh
		mydirectionalconvexfunction.hh
		quadraticenergy.hh
		tectonic.hh
		assemblers.hh
		explicitgrid.hh
		explicitvectors.hh
		gridselector.hh
		DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

src/assemblers.cc→dune/tectonic/assemblers.cc

+75 −64

Original line number	Diff line number	Diff line
		@@ -4,7 +4,8 @@

		#include <dune/istl/scaledidmatrix.hh>

		#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
		#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
		#include <dune/fufem/assemblers/boundaryoperatorassembler.hh>
		#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
		#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
		#include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
		@@ -12,41 +13,49 @@
		#include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
		#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
		#include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
		#include <dune/fufem/boundarypatch.hh>
		#include <dune/fufem/functions/basisgridfunction.hh>
		#include <dune/fufem/functions/constantfunction.hh>
		#include <dune/fufem/functiontools/p0p1interpolation.hh>
		#include <dune/fufem/quadraturerules/quadraturerulecache.hh>

		#include <dune/tectonic/frictionpotential.hh>
		#include <dune/tectonic/globalratestatefriction.hh>
		#include "data-structures/friction/frictionpotential.hh"
		#include "data-structures/friction/globalratestatefriction.hh"

		#include "assemblers.hh"

		template <class GridView, int dimension>
		MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
		: cellBasis(_gridView),
		MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView) :
		cellBasis(_gridView),
		vertexBasis(_gridView),
		gridView(_gridView),
		cellAssembler(cellBasis, cellBasis),
		vertexAssembler(vertexBasis, vertexBasis) {}

		template <class GridView, int dimension>
		template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
		void MyAssembler<GridView, dimension>::assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
		GlobalVectorType& b,
		const BoundaryPatch<GridView>& boundaryPatch,
		bool initializeVector) const {

		vertexAssembler.assembleBoundaryFunctional(localAssembler, b, boundaryPatch, initializeVector);
		}

		template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
		BoundaryPatch<GridView> const &frictionalBoundary,
		ScalarMatrix &frictionalBoundaryMass) const {
		BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
		LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
		frictionalBoundaryMassAssembler(frictionalBoundary);
		vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
		frictionalBoundaryMass);

		MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis> localMassAssembler;
		const BoundaryOperatorAssembler<VertexBasis, VertexBasis> boundaryAssembler(vertexBasis, vertexBasis, frictionalBoundary);
		boundaryAssembler.assemble(localMassAssembler, frictionalBoundaryMass);
		}

		template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleMass(
		Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
		densityFunction,
		Dune::VirtualFunction<LocalVector, LocalScalarVector> const & densityFunction,
		Matrix &M) const {

		// NOTE: We treat the weight as a constant function
		QuadratureRuleKey quadKey(dimension, 0);

		@@ -58,8 +67,11 @@ void MyAssembler<GridView, dimension>::assembleMass(
		}

		template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
		void MyAssembler<GridView, dimension>::assembleElasticity(
		double E,
		double nu,
		Matrix &A) const {

		StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
		localStiffness(E, nu);
		vertexAssembler.assembleOperator(localStiffness, A);
		@@ -70,6 +82,7 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
		Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
		Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
		Matrix &C) const {

		// NOTE: We treat the weights as constant functions
		QuadratureRuleKey shearViscosityKey(dimension, 0);
		QuadratureRuleKey bulkViscosityKey(dimension, 0);
		@@ -85,6 +98,7 @@ template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleBodyForce(
		Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
		Vector &f) const {

		L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
		gravityFunctionalAssembler(gravityField);
		vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
		@@ -92,11 +106,33 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(

		template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleNeumann(
		BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
		Dune::VirtualFunction<double, double> const &neumann,
		const BoundaryPatch<GridView>& neumannBoundary,
		const Dune::BitSetVector<dimension>& neumannNodes,
		Vector& f,
		const Dune::VirtualFunction<double, double>& neumann,
		double relativeTime) const {

		typename LocalVector::field_type val = 0;
		neumann.evaluate(relativeTime, val);

		size_t idx = 0;
		bool neumannIdx = false;
		for (; idx<neumannNodes.size() && !neumannIdx; idx++) {
		for (size_t d=0; d<dimension; d++) {
		if (neumannNodes[idx][d]) {
		neumannIdx = true;
		break;
		}
		}
		}
		idx--;

		LocalVector localNeumann(0);
		neumann.evaluate(relativeTime, localNeumann[0]);
		for (size_t i=0; i<localNeumann.size(); i++) {
		if (neumannNodes[idx][i])
		localNeumann[i] = val;
		}

		NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
		std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
		localNeumann));
		@@ -107,8 +143,11 @@ void MyAssembler<GridView, dimension>::assembleNeumann(
		template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
		BoundaryPatch<GridView> const &frictionalBoundary,
		ScalarVector &weightedNormalStress, double youngModulus,
		double poissonRatio, Vector const &displacement) const {
		ScalarVector &weightedNormalStress,
		ScalarVector &weights,
		double youngModulus,
		double poissonRatio,
		Vector const &displacement) const {

		BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
		displacement);
		@@ -121,7 +160,6 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
		auto const nodalTractionAverage =
		interpolateP0ToP1(frictionalBoundary, traction);

		ScalarVector weights;
		{
		NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
		frictionalBoundaryAssembler(
		@@ -130,47 +168,20 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
		vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
		weights, frictionalBoundary);
		}

		auto const normals = frictionalBoundary.getNormals();
		for (size_t i = 0; i < vertexBasis.size(); ++i)
		weightedNormalStress[i] =
		std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
		}

		template <class GridView, int dimension>
		auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
		Config::FrictionModel frictionModel,
		BoundaryPatch<GridView> const &frictionalBoundary,
		GlobalFrictionData<dimension> const &frictionInfo,
		ScalarVector const &weightedNormalStress) const
		-> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
		// Lumping of the nonlinearity
		ScalarVector weights;
		{
		NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
		frictionalBoundaryAssembler(std::make_shared<
		ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
		1));
		vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
		weights, frictionalBoundary);
		}
		switch (frictionModel) {
		case Config::Truncated:
		return std::make_shared<GlobalRateStateFriction<
		Matrix, Vector, TruncatedRateState, GridView>>(
		frictionalBoundary, frictionInfo, weights, weightedNormalStress);
		case Config::Regularised:
		return std::make_shared<GlobalRateStateFriction<
		Matrix, Vector, RegularisedRateState, GridView>>(
		frictionalBoundary, frictionInfo, weights, weightedNormalStress);
		default:
		assert(false);
		}
		}

		template <class GridView, int dimension>
		void MyAssembler<GridView, dimension>::assembleVonMisesStress(
		double youngModulus, double poissonRatio, Vector const &u,
		double youngModulus,
		double poissonRatio,
		Vector const &u,
		ScalarVector &stress) const {

		auto const gridDisplacement =
		std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
		vertexBasis, u);

dune/tectonic/assemblers.hh

0 → 100644

+110 −0

Original line number	Diff line number	Diff line
		#ifndef SRC_ASSEMBLERS_HH
		#define SRC_ASSEMBLERS_HH

		#include <dune/common/bitsetvector.hh>
		#include <dune/common/function.hh>
		#include <dune/istl/bcrsmatrix.hh>
		#include <dune/istl/bvector.hh>

		#include <dune/fufem/assemblers/assembler.hh>
		#pragma clang diagnostic push
		#pragma clang diagnostic ignored "-Wsign-compare"
		#include <dune/fufem/functionspacebases/p0basis.hh>
		#pragma clang diagnostic pop
		#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
		#include <dune/fufem/boundarypatch.hh>

		#include "data-structures/friction/globalfriction.hh"
		#include "data-structures/friction/globalfrictiondata.hh"

		#include "data-structures/enums.hh"

		template <class GridView, int dimension>
		class MyAssembler {
		public:
		using GV = GridView;
		using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
		using Matrix =
		Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
		using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
		using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;

		using CellBasis = P0Basis<GridView, double>;
		using VertexBasis = P1NodalBasis<GridView, double>;

		CellBasis const cellBasis;
		VertexBasis const vertexBasis;
		GridView const &gridView;

		private:
		using Grid = typename GridView::Grid;
		using LocalVector = typename Vector::block_type;
		using LocalScalarVector = typename ScalarVector::block_type;

		using LocalCellBasis = typename CellBasis::LocalFiniteElement;
		using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;

		Assembler<CellBasis, CellBasis> cellAssembler;
		Assembler<VertexBasis, VertexBasis> vertexAssembler;

		public:
		MyAssembler(GridView const &gridView);

		template <class LocalBoundaryFunctionalAssemblerType, class GlobalVectorType>
		void assembleBoundaryFunctional(LocalBoundaryFunctionalAssemblerType& localAssembler,
		GlobalVectorType& b,
		const BoundaryPatch<GridView>& boundaryPatch,
		bool initializeVector=true) const;

		void assembleFrictionalBoundaryMass(
		BoundaryPatch<GridView> const &frictionalBoundary,
		ScalarMatrix &frictionalBoundaryMass) const;

		void assembleMass(
		Dune::VirtualFunction<
		LocalVector, LocalScalarVector> const &densityFunction,
		Matrix &M) const;

		void assembleElasticity(
		double E,
		double nu,
		Matrix &A) const;

		void assembleViscosity(
		Dune::VirtualFunction<LocalVector, LocalScalarVector> const & shearViscosity,
		Dune::VirtualFunction<LocalVector, LocalScalarVector> const & bulkViscosity,
		Matrix &C) const;

		void assembleBodyForce(
		Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
		Vector &f) const;

		void assembleNeumann(
		const BoundaryPatch<GridView>& neumannBoundary,
		const Dune::BitSetVector<dimension>& neumannNodes,
		Vector& f,
		const Dune::VirtualFunction<double, double>& neumann,
		double relativeTime) const;

		void assembleWeightedNormalStress(
		BoundaryPatch<GridView> const &frictionalBoundary,
		ScalarVector &weightedNormalStress,
		ScalarVector &weights,
		double youngModulus,
		double poissonRatio,
		Vector const &displacement) const;

		auto assembleFrictionNonlinearity(
		Config::FrictionModel frictionModel,
		BoundaryPatch<GridView> const &frictionalBoundary,
		GlobalFrictionData<dimension> const &frictionInfo,
		ScalarVector const &weightedNormalStress) const
		-> std::shared_ptr<GlobalFriction<Matrix, Vector>>;

		void assembleVonMisesStress(
		double youngModulus,
		double poissonRatio,
		Vector const &u,
		ScalarVector &stress) const;
		};
		#endif

dune/tectonic/assemblers_tmpl.cc

0 → 100644

+22 −0

Original line number	Diff line number	Diff line
		#ifndef MY_DIM
		#error MY_DIM unset
		#endif

		#include "explicitgrid.hh"
		#include "explicitvectors.hh"

		#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
		#include <dune/fufem/boundarypatch.hh>

		#include "assemblers.hh"

		template class MyAssembler<DefLeafGridView, MY_DIM>;


		using MyNeumannBoundaryAssembler = NeumannBoundaryAssembler<DeformedGrid, typename ScalarVector::block_type>;

		template void MyAssembler<DefLeafGridView, MY_DIM>::assembleBoundaryFunctional<MyNeumannBoundaryAssembler, ScalarVector>(
		MyNeumannBoundaryAssembler& localAssembler,
		ScalarVector& b,
		const BoundaryPatch<DefLeafGridView>& boundaryPatch,
		bool initializeVector) const;

dune/tectonic/data-structures/.cc

0 → 100644

+7 −0

Original line number	Diff line number	Diff line
		#ifndef MY_DIM
		#error MY_DIM unset
		#endif

		#include "../explicitgrid.hh"

		template class ContactNetwork<Grid, MY_DIM>;

dune/tectonic/data-structures/CMakeLists.txt

0 → 100644

+19 −0

Original line number	Diff line number	Diff line
		add_subdirectory("body")
		add_subdirectory("friction")
		add_subdirectory("network")

		add_custom_target(tectonic_dune_data-structures SOURCES
		enumparser.hh
		enumparser.cc
		enums.hh
		matrices.hh
		program_state.hh
		)

		#install headers
		install(FILES
		enumparser.hh
		enums.hh
		matrices.hh
		program_state.hh
		DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/data-structures/body/CMakeLists.txt

0 → 100644

+13 −0

Original line number	Diff line number	Diff line
		add_custom_target(tectonic_dune_data-structures_body SOURCES
		body.hh
		body.cc
		bodydata.hh
		boundarycondition.hh
		)

		#install headers
		install(FILES
		body.hh
		bodydata.hh
		boundarycondition.hh
		DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/data-structures/body/body.cc

0 → 100644

+347 −0

Original line number	Diff line number	Diff line
		#ifdef HAVE_CONFIG_H
		#include "config.h"
		#endif

		#include <dune/common/hybridutilities.hh>
		#include "body.hh"

		#include <dune/tectonic/utils/debugutils.hh>

		// -------------------------------
		// --- LeafBody Implementation ---
		// -------------------------------
		template <class GridTEMPLATE, class VectorType>
		LeafBody<GridTEMPLATE, VectorType>::LeafBody(
		const std::shared_ptr<BodyData<dim>>& bodyData,
		const std::shared_ptr<GridTEMPLATE>& hostGrid) :
		bodyData_(bodyData),
		hostGrid_(hostGrid),
		deformation_(std::make_shared<DeformationFunction>(hostGrid_->leafGridView())),
		grid_(std::make_shared<GridType>(hostGrid_, deformation_)),
		gridView_(grid_->leafGridView()),
		assembler_(std::make_shared<Assembler>(gridView_)),
		matrices_()
		{

		boundaryConditions_.resize(0);

		externalForce_ = [&](double relativeTime, VectorType &ell) {
		// Problem formulation: right-hand side
		std::vector<std::shared_ptr<BoundaryCondition>> leafNeumannConditions;
		boundaryConditions("neumann", leafNeumannConditions);
		ell.resize(gravityFunctional_.size());
		ell = gravityFunctional_;

		for (size_t i=0; i<leafNeumannConditions.size(); i++) {
		const auto& leafNeumannCondition = leafNeumannConditions[i];

		VectorType b;
		assembler_->assembleNeumann(*leafNeumannCondition->boundaryPatch(),
		*leafNeumannCondition->boundaryNodes(),
		b,
		*leafNeumannCondition->boundaryFunction(),
		relativeTime);
		ell += b;
		}
		};
		}


		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::assemble() {

		// assemble matrices
		assembler_->assembleElasticity(bodyData_->getYoungModulus(), bodyData_->getPoissonRatio(), *matrices_.elasticity);
		//print(*matrices_.elasticity, "elasticity");

		assembler_->assembleViscosity(bodyData_->getShearViscosityField(), bodyData_->getBulkViscosityField(), *matrices_.damping);
		//print(*matrices_.damping, "viscosity");

		assembler_->assembleMass(bodyData_->getDensityField(), *matrices_.mass);
		//print(*matrices_.mass, "mass");

		// assemble forces
		assembler_->assembleBodyForce(bodyData_->getGravityField(), gravityFunctional_);
		//print(gravityFunctional_, "gravity");
		}

		// setter and getter

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::data() const
		-> const std::shared_ptr<BodyData<dim>>& {

		return bodyData_;
		}


		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::setDeformation(const VectorType& def) {
		deformation_->setDeformation(def);
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::deformation() const
		-> DeformationFunction& {

		return *deformation_;
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::grid() const
		-> std::shared_ptr<GridType> {

		return grid_;
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::hostGrid() const
		-> std::shared_ptr<GridTEMPLATE> {

		return hostGrid_;
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::gridView() const
		-> const GridView& {

		return gridView_;
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::nVertices() const
		-> size_t {
		return gridView_.size(dim);
		}


		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::assembler() const
		-> const std::shared_ptr<Assembler>& {

		return assembler_;
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::matrices() const
		-> const Matrices<Matrix,1>& {

		return matrices_;
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::externalForce() const
		-> const ExternalForce& {

		return externalForce_;
		}

		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {
		boundaryConditions_.emplace_back(boundaryCondition);
		}

		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::boundaryConditions(
		const std::string& tag,
		std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {

		selectedConditions.resize(0);
		for (size_t i=0; i<boundaryConditions_.size(); i++) {
		if (boundaryConditions_[i]->tag() == tag)
		selectedConditions.emplace_back(boundaryConditions_[i]);
		}
		}

		template <class GridTEMPLATE, class VectorType>
		auto LeafBody<GridTEMPLATE, VectorType>::boundaryConditions() const
		-> const std::vector<std::shared_ptr<BoundaryCondition>>& {

		return boundaryConditions_;
		}

		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::boundaryPatchNodes(
		const std::string& tag,
		BoundaryPatchNodes& nodes) const {

		std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		nodes.resize(_boundaryConditions.size());

		for (size_t i=0; i<nodes.size(); i++) {
		nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
		}
		}

		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::boundaryNodes(
		const std::string& tag,
		BoundaryNodes& nodes) const {

		std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		nodes.resize(_boundaryConditions.size());

		for (size_t i=0; i<nodes.size(); i++) {
		nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
		}
		}

		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::boundaryFunctions(
		const std::string& tag,
		BoundaryFunctions& functions) const {

		std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		functions.resize(_boundaryConditions.size());

		for (size_t i=0; i<functions.size(); i++) {
		functions[i] = _boundaryConditions[i]->boundaryFunction().get();
		}
		}

		template <class GridTEMPLATE, class VectorType>
		void LeafBody<GridTEMPLATE, VectorType>::boundaryPatches(
		const std::string& tag,
		BoundaryPatches& patches) const {

		std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		patches.resize(_boundaryConditions.size());

		for (size_t i=0; i<patches.size(); i++) {
		patches[i] = _boundaryConditions[i]->boundaryPatch().get();
		}
		}


		// ---------------------------
		// --- Body Implementation ---
		// ---------------------------

		// setter and getter

		template <class GridView>
		auto Body<GridView>::data() const
		-> const std::shared_ptr<BodyData<dim>>& {

		return bodyData_;
		}


		template <class GridView>
		auto Body<GridView>::grid() const
		-> std::shared_ptr<GridType> {

		return grid_;
		}

		template <class GridView>
		auto Body<GridView>::gridView() const
		-> const GridView& {

		return gridView_;
		}

		template <class GridView>
		auto Body<GridView>::nVertices() const
		-> size_t {
		return gridView_.size(dim);
		}


		template <class GridView>
		void Body<GridView>::addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition) {

		boundaryConditions_.push_back(boundaryCondition);
		}


		template <class GridView>
		void Body<GridView>::boundaryConditions(
		const std::string& tag,
		std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const {

		selectedConditions.resize(0);
		for (size_t i=0; i<boundaryConditions_.size(); i++) {
		if (boundaryConditions_[i]->tag() == tag)
		selectedConditions.push_back(boundaryConditions_[i]);
		}
		}


		template <class GridView>
		auto Body<GridView>::boundaryConditions() const
		-> const std::vector<std::shared_ptr<BoundaryCondition>>& {

		return boundaryConditions_;
		}


		template <class GridView>
		void Body<GridView>::boundaryPatchNodes(
		const std::string& tag,
		BoundaryPatchNodes& nodes) const {

		std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		nodes.resize(_boundaryConditions.size());

		for (size_t i=0; i<nodes.size(); i++) {
		nodes[i] = _boundaryConditions[i]->boundaryPatch()->getVertices();
		}
		}

		template <class GridView>
		void Body<GridView>::boundaryNodes(
		const std::string& tag,
		BoundaryNodes& nodes) const {

		std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		nodes.resize(_boundaryConditions.size());

		for (size_t i=0; i<nodes.size(); i++) {
		nodes[i] = _boundaryConditions[i]->boundaryNodes().get();
		}
		}

		template <class GridView>
		void Body<GridView>::boundaryFunctions(
		const std::string& tag,
		BoundaryFunctions& functions) const {

		std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		functions.resize(_boundaryConditions.size());

		for (size_t i=0; i<functions.size(); i++) {
		functions[i] = _boundaryConditions[i]->boundaryFunction().get();
		}
		}

		template <class GridView>
		void Body<GridView>::boundaryPatches(
		const std::string& tag,
		BoundaryPatches& patches) const {

		std::vector<std::shared_ptr<BoundaryCondition>> _boundaryConditions;
		this->boundaryConditions(tag, _boundaryConditions);

		patches.resize(_boundaryConditions.size());

		for (size_t i=0; i<patches.size(); i++) {
		patches[i] = _boundaryConditions[i]->boundaryPatch().get();
		}
		}

		#include "body_tmpl.cc"

dune/tectonic/data-structures/body/body.hh

0 → 100644

+209 −0

Original line number	Diff line number	Diff line
		#ifndef SRC_DATA_STRUCTURES_BODY_HH
		#define SRC_DATA_STRUCTURES_BODY_HH

		#include <functional>
		#include <type_traits>

		#include <dune/common/bitsetvector.hh>
		#include <dune/common/function.hh>
		#include <dune/common/shared_ptr.hh>

		#include <dune/istl/bcrsmatrix.hh>
		#include <dune/istl/bvector.hh>

		//#include <dune/fufem/assemblers/assembler.hh>
		//#pragma clang diagnostic push
		//#pragma clang diagnostic ignored "-Wsign-compare"
		//#include <dune/fufem/functionspacebases/p0basis.hh>
		//#pragma clang diagnostic pop
		//#include <dune/fufem/functionspacebases/p1nodalbasis.hh>

		#include <dune/grid/geometrygrid/grid.hh>

		#include <dune/fufem/functions/deformationfunction.hh>

		#include <dune/solvers/norms/energynorm.hh>

		#include "../../assemblers.hh"
		#include "../enums.hh"
		#include "../matrices.hh"

		#include "boundarycondition.hh"
		#include "bodydata.hh"

		template <class HostGridType, class VectorType>
		class LeafBody {
		public:
		static const int dim = HostGridType::dimensionworld;

		using DeformationFunction = DeformationFunction<typename HostGridType::LeafGridView, VectorType>;
		using GridType = Dune::GeometryGrid<HostGridType, DeformationFunction>;
		using GridView = typename GridType::LeafGridView;

		using Function = Dune::VirtualFunction<double, double>;
		using ExternalForce = std::function<void(double, VectorType &)>;

		using Assembler = MyAssembler<GridView, dim>;
		using Matrix = typename Assembler::Matrix;
		using LocalVector = typename VectorType::block_type;

		using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
		using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;

		using BoundaryCondition = BoundaryCondition<GridView, dim>;

		using BoundaryFunctions = std::vector<const Function* >;
		using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
		using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
		using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;

		using StateEnergyNorm = EnergyNorm<ScalarMatrix, ScalarVector>;

		private:
		std::shared_ptr<BodyData<dim>> bodyData_;

		std::shared_ptr<HostGridType> hostGrid_;
		std::shared_ptr<DeformationFunction> deformation_;
		std::shared_ptr<GridType> grid_;
		GridView gridView_;

		std::shared_ptr<Assembler> assembler_;
		Matrices<Matrix,1> matrices_;

		ExternalForce externalForce_;

		VectorType gravityFunctional_;

		// boundary conditions
		std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;

		public:
		LeafBody(
		const std::shared_ptr<BodyData<dim>>& bodyData,
		const std::shared_ptr<HostGridType>& hostGrid);

		void assemble();

		// setter and getter
		auto data() const -> const std::shared_ptr<BodyData<dim>>&;

		void setDeformation(const VectorType& def);
		auto deformation() const -> DeformationFunction&;
		auto grid() const -> std::shared_ptr<GridType>;
		auto hostGrid() const -> std::shared_ptr<HostGridType>;

		auto gridView() const -> const GridView&;
		auto nVertices() const -> size_t;

		auto assembler() const -> const std::shared_ptr<Assembler>&;
		auto matrices() const -> const Matrices<Matrix,1>&;

		auto externalForce() const -> const ExternalForce&;

		void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
		void boundaryConditions(
		const std::string& tag,
		std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
		auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;

		void boundaryPatchNodes(
		const std::string& tag,
		BoundaryPatchNodes& nodes) const;
		void boundaryNodes(
		const std::string& tag,
		BoundaryNodes& nodes) const;
		void boundaryFunctions(
		const std::string& tag,
		BoundaryFunctions& functions) const;
		void boundaryPatches(
		const std::string& tag,
		BoundaryPatches& patches) const;
		};

		template <class GridView>
		class Body {
		public:
		enum {dim = GridView::dimensionworld};

		using GridType = typename GridView::Grid;

		using Function = Dune::VirtualFunction<double, double>;

		using BoundaryCondition = BoundaryCondition<GridView, dim>;

		using BoundaryFunctions = std::vector<const Function* >;
		using BoundaryNodes = std::vector<const Dune::BitSetVector<dim>* >;
		using BoundaryPatchNodes = std::vector<const Dune::BitSetVector<1>* >;
		using BoundaryPatches = std::vector<const typename BoundaryCondition::BoundaryPatch* >;

		private:
		std::shared_ptr<BodyData<dim>> bodyData_;

		std::shared_ptr<GridType> grid_;
		const size_t level_;
		GridView gridView_;

		// boundary conditions
		std::vector<std::shared_ptr<BoundaryCondition>> boundaryConditions_;

		public:
		template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LeafGridView>::value, int> = 0>
		Body(
		const std::shared_ptr<BodyData<dim>>& bodyData,
		const std::shared_ptr<GridType>& grid) :
		bodyData_(bodyData),
		grid_(grid),
		level_(grid_->maxLevel()),
		gridView_(grid_->leafGridView()) {

		boundaryConditions_.resize(0);
		}

		template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
		Body(
		const std::shared_ptr<BodyData<dim>>& bodyData,
		const std::shared_ptr<GridType>& grid,
		const size_t level) :
		bodyData_(bodyData),
		grid_(grid),
		level_(level),
		gridView_(grid_->levelGridView(level_)) {

		boundaryConditions_.resize(0);
		}


		// setter and getter
		auto data() const -> const std::shared_ptr<BodyData<dim>>&;
		auto grid() const -> std::shared_ptr<GridType>;

		//template <class GridView_ = GridView, std::enable_if_t<std::is_same<GridView_, typename GridType::LevelGridView>::value, int> = 0>
		auto level() const
		-> size_t {

		return level_;
		}

		auto gridView() const -> const GridView&;
		auto nVertices() const -> size_t;

		void addBoundaryCondition(std::shared_ptr<BoundaryCondition> boundaryCondition);
		void boundaryConditions(
		const std::string& tag,
		std::vector<std::shared_ptr<BoundaryCondition>>& selectedConditions) const;
		auto boundaryConditions() const -> const std::vector<std::shared_ptr<BoundaryCondition>>&;

		void boundaryPatchNodes(
		const std::string& tag,
		BoundaryPatchNodes& nodes) const;
		void boundaryNodes(
		const std::string& tag,
		BoundaryNodes& nodes) const;
		void boundaryFunctions(
		const std::string& tag,
		BoundaryFunctions& functions) const;
		void boundaryPatches(
		const std::string& tag,
		BoundaryPatches& patches) const;
		};
		#endif

dune/tectonic/data-structures/body/body_tmpl.cc

0 → 100644

+11 −0

Original line number	Diff line number	Diff line
		#ifndef MY_DIM
		#error MY_DIM unset
		#endif

		#include "../../explicitgrid.hh"
		#include "../../explicitvectors.hh"

		template class LeafBody<Grid, Vector>;

		template class Body<DefLeafGridView>;
		template class Body<DefLevelGridView>;

dune/tectonic/body.hh→dune/tectonic/data-structures/body/bodydata.hh

+3 −3

Original line number	Diff line number	Diff line
		#ifndef DUNE_TECTONIC_BODY_HH
		#define DUNE_TECTONIC_BODY_HH
		#ifndef DUNE_TECTONIC_BODY_DATA_HH
		#define DUNE_TECTONIC_BODY_DATA_HH

		template <int dimension> struct Body {
		template <int dimension> struct BodyData {
		using ScalarFunction =
		Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
		Dune::FieldVector<double, 1>>;

dune/tectonic/data-structures/body/boundarycondition.hh

0 → 100644

+87 −0

Original line number	Diff line number	Diff line
		#ifndef SRC_BOUNDARYCONDITION_HH
		#define SRC_BOUNDARYCONDITION_HH

		#include <dune/common/bitsetvector.hh>
		#include <dune/common/function.hh>

		#include <dune/fufem/boundarypatch.hh>

		#include <dune/tectonic/utils/tobool.hh>

		template <class GridView, int dims>
		class BoundaryCondition {
		public:
		using BoundaryPatch = BoundaryPatch<GridView>;

		private:
		using Function = Dune::VirtualFunction<double, double>;

		const std::string tag_; // store type of boundary condition, e.g. dirichlet, neumann, friction, etc

		std::shared_ptr<BoundaryPatch> boundaryPatch_;
		std::shared_ptr<Dune::BitSetVector<dims>> boundaryNodes_;
		std::shared_ptr<Function> boundaryFunction_;

		public:
		BoundaryCondition(const std::string& tag = "") :
		tag_(tag)
		{}

		BoundaryCondition(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function, const std::string& tag = "") :
		tag_(tag),
		boundaryFunction_(function)
		{
		setBoundaryPatch(patch);
		}

		void setBoundaryPatch(const GridView& gridView, std::shared_ptr<Dune::BitSetVector<dims>> nodes) {
		boundaryNodes_ = nodes;
		boundaryPatch_ = std::make_shared<BoundaryPatch>(gridView, *nodes);
		}

		void setBoundaryPatch(std::shared_ptr<BoundaryPatch> patch) {
		boundaryPatch_ = patch;

		auto nodes = patch->getVertices();
		boundaryNodes_ = std::make_shared<Dune::BitSetVector<dims>>(nodes->size(), false);
		for (size_t i=0; i<nodes->size(); i++) {
		if (toBool((*nodes)[i])) {
		for (size_t d=0; d<dims; d++) {
		(*boundaryNodes_)[i][d] = true;
		}
		}
		}
		}

		void setBoundaryPatch(const BoundaryPatch& patch) {
		auto patch_ptr = std::make_shared<BoundaryPatch>(patch);
		setBoundaryPatch(patch_ptr);
		}

		void setBoundaryFunction(std::shared_ptr<Function> function) {
		boundaryFunction_ = function;
		}

		void set(std::shared_ptr<BoundaryPatch> patch, std::shared_ptr<Function> function) {
		setBoundaryPatch(patch);
		boundaryFunction_ = function;
		}

		const std::string& tag() const {
		return tag_;
		}

		const std::shared_ptr<BoundaryPatch>& boundaryPatch() const {
		return boundaryPatch_;
		}

		const std::shared_ptr<Dune::BitSetVector<dims>>& boundaryNodes() const {
		return boundaryNodes_;
		}

		const std::shared_ptr<Function>& boundaryFunction() const {
		return boundaryFunction_;
		}
		};

		#endif

src/enumparser.cc→dune/tectonic/data-structures/enumparser.cc

+6 −0

Original line number	Diff line number	Diff line
		@@ -32,6 +32,12 @@ Config::FrictionModel StringToEnum<Config::FrictionModel>::convert(
		if (s == "Regularised")
		return Config::Regularised;

		if (s == "Tresca")
		return Config::Tresca;

		if (s == "None")
		return Config::None;

		DUNE_THROW(Dune::Exception, "failed to parse enum");
		}

src/enumparser.hh→dune/tectonic/data-structures/enumparser.hh

+2 −2

Original line number	Diff line number	Diff line
		#ifndef SRC_ENUMPARSER_HH
		#define SRC_ENUMPARSER_HH
		#ifndef SRC_DATA_STRUCTURES_ENUMPARSER_HH
		#define SRC_DATA_STRUCTURES_ENUMPARSER_HH

		// Copyright Carsten Graeser 2012

src/enums.hh→dune/tectonic/data-structures/enums.hh

+1 −1

Original line number	Diff line number	Diff line
		@@ -2,7 +2,7 @@
		#define SRC_ENUMS_HH

		struct Config {
		enum FrictionModel { Truncated, Regularised };
		enum FrictionModel { Truncated, Regularised, Tresca, None };
		enum stateModel { AgeingLaw, SlipLaw };
		enum scheme { Newmark, BackwardEuler };
		enum PatchType { Rectangular, Trapezoidal };

dune/tectonic/data-structures/friction/CMakeLists.txt

0 → 100644

+20 −0

Original line number	Diff line number	Diff line
		add_custom_target(tectonic_dune_data-structures_friction SOURCES
		frictioncouplingpair.hh
		frictiondata.hh
		frictionpotential.hh
		globalfriction.hh
		globalfrictiondata.hh
		globalratestatefriction.hh
		localfriction.hh
		)

		#install headers
		install(FILES
		frictioncouplingpair.hh
		frictiondata.hh
		frictionpotential.hh
		globalfriction.hh
		globalfrictiondata.hh
		globalratestatefriction.hh
		localfriction.hh
		DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/data-structures/friction/frictioncouplingpair.hh

0 → 100644

+38 −0

Original line number	Diff line number	Diff line
		#ifndef SRC_FRICTIONCOUPLINGPAIR_HH
		#define SRC_FRICTIONCOUPLINGPAIR_HH

		#include <dune/fufem/geometry/convexpolyhedron.hh>
		#include <dune/contact/common/couplingpair.hh>

		#include "globalfrictiondata.hh"

		template <class GridType, class LocalVectorType, class field_type = double>
		class FrictionCouplingPair : public Dune::Contact::CouplingPair<GridType, GridType, field_type>{
		private:
		static const int dim = GridType::dimensionworld;

		using Base = Dune::Contact::CouplingPair<GridType,GridType,field_type>;
		using LocalVector = LocalVectorType;

		// friction data
		std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch_;
		std::shared_ptr<GlobalFrictionData<dim>> frictionData_;

		public:
		void setWeakeningPatch(std::shared_ptr<ConvexPolyhedron<LocalVector>> weakeningPatch) {
		weakeningPatch_ = weakeningPatch;
		}

		void setFrictionData(std::shared_ptr<GlobalFrictionData<dim>> frictionData) {
		frictionData_ = frictionData;
		}

		const auto& weakeningPatch() const {
		return *weakeningPatch_;
		}

		const auto& frictionData() const {
		return *frictionData_;
		}
		};
		#endif

dune/tectonic/frictiondata.hh→dune/tectonic/data-structures/friction/frictiondata.hh

+0 −0

File moved.

dune/tectonic/frictionpotential.hh→dune/tectonic/data-structures/friction/frictionpotential.hh

+72 −5

Original line number	Diff line number	Diff line
		@@ -9,7 +9,7 @@
		#include <dune/common/exceptions.hh>
		#include <dune/common/function.hh>

		#include <dune/tectonic/frictiondata.hh>
		#include "frictiondata.hh"

		class FrictionPotential {
		public:
		@@ -34,33 +34,56 @@ class TruncatedRateState : public FrictionPotential {
		: fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}

		double coefficientOfFriction(double V) const override {
		if (V <= Vmin)
		if (V <= Vmin or regularity(V)>10e8)
		return 0.0;

		return fd.a * std::log(V / Vmin);
		//std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;

		auto res = fd.a * std::log(V / Vmin);
		if (std::isinf(res)) {
		//std::cout << "V: " << V << " Vmin: " << Vmin << std::endl;
		}

		return res;
		}

		double differential(double V) const override {
		//std::cout << "differential: " << weight * fd.C - weightedNormalStress * coefficientOfFriction(V) << std::endl;
		if (V <= Vmin or regularity(V)>10e8)
		return 0.0;

		return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
		}

		double second_deriv(double V) const override {
		//std::cout << "second deriv: " << ((V<=Vmin) ? 0.0 : -weightedNormalStress * (fd.a / V)) << std::endl;

		if (V <= Vmin)
		return 0;
		return 0.0;

		return - weightedNormalStress * (fd.a / V);
		}

		double regularity(double V) const override {
		//std::cout << "regularity: " << ((std::abs(V - Vmin) < 1e-14) ? std::numeric_limits<double>::infinity() : std::abs(second_deriv(V))) << std::endl;

		if (std::abs(V - Vmin) < 1e-14) // TODO
		return std::numeric_limits<double>::infinity();

		return std::abs(second_deriv(V));
		}

		double evaluate(double V) const override {
		if (V <= Vmin or regularity(V)>10e8)
		return 0.0;

		return weight * fd.C * V - weightedNormalStress * fd.a * V * (std::log(V / Vmin) - 1);
		}

		void updateAlpha(double alpha) override {
		double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
		Vmin = fd.V0 / std::exp(logrest);
		//Vmin = 0.0;
		}

		private:
		@@ -92,6 +115,10 @@ class RegularisedRateState : public FrictionPotential {
		return std::abs(second_deriv(V));
		}

		double evaluate(double V) const override {
		return weight * fd.C * V - weightedNormalStress * 4 * fd.a * Vmin * std::pow(std::asinh(0.25 * V / Vmin), 2.0);
		}

		void updateAlpha(double alpha) override {
		double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
		Vmin = fd.V0 / std::exp(logrest);
		@@ -104,8 +131,48 @@ class RegularisedRateState : public FrictionPotential {
		double Vmin;
		};

		class Tresca : public FrictionPotential {
		public:
		Tresca(double _weight, double _weightedNormalStress, FrictionData _fd)
		: fd(_fd), weightedNormalStress(_weightedNormalStress) {}

		double coefficientOfFriction(double V) const override {
		return fd.mu0;
		}

		double differential(double V) const override {
		//if (std::isinf(regularity(V)))
		// return 0.0;
		return - weightedNormalStress * coefficientOfFriction(V);
		}

		double second_deriv(double V) const override {
		return 0;
		}

		double regularity(double V) const override {
		if (std::abs(V) < 1e-14) // TODO
		return std::numeric_limits<double>::infinity();

		return std::abs(second_deriv(V));
		}

		double evaluate(double V) const override {
		return - weightedNormalStress * fd.mu0 * V;
		}

		void updateAlpha(double alpha) override {}

		private:
		FrictionData const fd;
		double const weightedNormalStress;
		};

		class ZeroFunction : public FrictionPotential {
		public:
		template <typename... Args>
		ZeroFunction(Args... args) {}

		double evaluate(double) const override { return 0; }

		double coefficientOfFriction(double s) const override { return 0; }

dune/tectonic/globalfriction.hh→dune/tectonic/data-structures/friction/globalfriction.hh

+9 −5

Original line number	Diff line number	Diff line
		@@ -9,9 +9,10 @@

		#include <dune/solvers/common/interval.hh>

		#include <dune/tectonic/localfriction.hh>
		#include "localfriction.hh"

		template <class Matrix, class Vector> class GlobalFriction {
		template <class Matrix, class Vector>
		class GlobalFriction {
		protected:
		using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;

		@@ -38,9 +39,10 @@ template <class Matrix, class Vector> class GlobalFriction {
		LocalNonlinearity const virtual &restriction(size_t i) const = 0;

		void addHessian(Vector const &v, Matrix &hessian) const {
		for (size_t i = 0; i < v.size(); ++i)
		for (size_t i = 0; i < v.size(); ++i) {
		restriction(i).addHessian(v[i], hessian[i][i]);
		}
		}

		void directionalDomain(Vector const &, Vector const &,
		Dune::Solvers::Interval<double> &dom) const {
		@@ -76,11 +78,13 @@ template <class Matrix, class Vector> class GlobalFriction {

		ScalarVector coefficientOfFriction(Vector const &x) const {
		ScalarVector ret(x.size());
		for (size_t i = 0; i < x.size(); ++i)
		for (size_t i = 0; i < x.size(); ++i) {
		ret[i] = restriction(i).coefficientOfFriction(x[i]);
		}
		return ret;
		}

		void virtual updateAlpha(ScalarVector const &alpha) = 0;
		void virtual updateAlpha(const std::vector<ScalarVector>& alpha) = 0;

		};
		#endif

dune/tectonic/globalfrictiondata.hh→dune/tectonic/data-structures/friction/globalfrictiondata.hh

+2 −1

Original line number	Diff line number	Diff line
		@@ -4,7 +4,7 @@
		#include <dune/common/function.hh>
		#include <dune/common/fvector.hh>

		#include <dune/tectonic/frictiondata.hh>
		#include "frictiondata.hh"

		template <class DomainType>
		double evaluateScalarFunction(
		@@ -27,6 +27,7 @@ template <int dimension> class GlobalFrictionData {
		Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
		Dune::FieldVector<double, 1>>;

		public:
		double virtual const &C() const = 0;
		double virtual const &L() const = 0;
		double virtual const &V0() const = 0;

dune/tectonic/data-structures/friction/globalratestatefriction.hh

0 → 100644

+127 −0

Original line number	Diff line number	Diff line
		#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
		#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH

		#include <vector>

		#include <dune/common/bitsetvector.hh>
		#include <dune/common/fmatrix.hh>
		#include <dune/common/fvector.hh>
		#include <dune/grid/common/mcmgmapper.hh>
		#include <dune/istl/bcrsmatrix.hh>
		#include <dune/istl/bvector.hh>

		#include "../../spatial-solving/contact/dualmortarcoupling.hh"

		#include "globalfrictiondata.hh"
		#include "globalfriction.hh"
		#include "frictioncouplingpair.hh"

		#include "../../utils/geocoordinate.hh"
		#include "../../utils/index-in-sorted-range.hh"

		template <class Matrix, class Vector, class ScalarFriction, class GridType>
		class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
		public:
		using GlobalFriction<Matrix, Vector>::block_size;
		using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;

		private:
		using Base = GlobalFriction<Matrix, Vector>;

		using field_type = typename Vector::field_type;
		using typename Base::ScalarVector;
		using typename Base::LocalVectorType;

		using FrictionCoupling = FrictionCouplingPair<GridType, LocalVectorType, field_type>;
		using ContactCoupling = DualMortarCoupling<field_type, GridType>;

		size_t bodyIndex(const size_t globalIdx) {
		size_t i = offSet_.size()-1;

		for (; i>0; ) {
		if (globalIdx >= offSet_[i])
		break;
		i--;
		}
		return i;
		}

		public:
		GlobalRateStateFriction(const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
		const std::vector<std::shared_ptr<FrictionCoupling>>& couplings, // contains frictionInfo
		const std::vector<ScalarVector>& weights,
		const std::vector<ScalarVector>& weightedNormalStress)
		: restrictions_(), localToGlobal_(), zeroFriction_() {

		assert(contactCouplings.size() == couplings.size());
		assert(weights.size() == weightedNormalStress.size());

		const auto nBodies = weights.size();
		offSet_.resize(nBodies, 0);
		for (size_t i=1; i<nBodies; i++) {
		offSet_[i] = offSet_[i-1] + weights[i-1].size();
		}

		std::vector<std::vector<int>> nonmortarBodies(nBodies); // first index body, second index coupling
		for (size_t i=0; i<contactCouplings.size(); i++) {
		const auto nonmortarIdx = couplings[i]->gridIdx_[0];
		nonmortarBodies[nonmortarIdx].emplace_back(i);
		}

		for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
		const auto& couplingIndices = nonmortarBodies[bodyIdx];

		if (couplingIndices.size()==0)
		continue;

		const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();

		Dune::MultipleCodimMultipleGeomTypeMapper<
		decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());

		for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
		const auto i = vertexMapper.index(*it);

		for (size_t j=0; j<couplingIndices.size(); j++) {
		const auto couplingIdx = couplingIndices[j];

		if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
		continue;

		localToGlobal_.emplace_back(offSet_[bodyIdx] + i);
		restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
		couplings[j]->frictionData()(geoToPoint(it->geometry())));
		break;
		}
		}
		}
		}

		void updateAlpha(const std::vector<ScalarVector>& alpha) override {
		//print(alpha, "alpha");
		for (size_t j = 0; j < restrictions_.size(); ++j) {
		const auto globalDof = localToGlobal_[j];
		const auto bodyIdx = bodyIndex(globalDof);
		size_t bodyDof = globalDof - offSet_[bodyIdx];

		restrictions_[j].updateAlpha(alpha[bodyIdx][bodyDof]);
		}
		}

		/*
		Return a restriction of the outer function to the i'th node.
		*/
		LocalNonlinearity const &restriction(size_t i) const override {
		auto const index = indexInSortedRange(localToGlobal_, i);
		if (index == localToGlobal_.size())
		return zeroFriction_;
		return restrictions_[index];
		}

		private:
		std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions_;
		std::vector<size_t> offSet_; // index off-set by body id
		std::vector<size_t> localToGlobal_;
		WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction_;
		};
		#endif

dune/tectonic/localfriction.hh→dune/tectonic/data-structures/friction/localfriction.hh

+66 −9

Original line number	Diff line number	Diff line
		@@ -10,7 +10,7 @@
		#include <dune/fufem/arithmetic.hh>
		#include <dune/solvers/common/interval.hh>

		#include <dune/tectonic/frictionpotential.hh>
		#include "frictionpotential.hh"

		template <size_t dimension> class LocalFriction {
		public:
		@@ -19,6 +19,7 @@ template <size_t dimension> class LocalFriction {
		using VectorType = Dune::FieldVector<double, dimension>;
		using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;

		double virtual operator()(VectorType const &x) const = 0;
		void virtual updateAlpha(double alpha) = 0;
		double virtual regularity(VectorType const &x) const = 0;
		double virtual coefficientOfFriction(VectorType const &x) const = 0;
		@@ -39,34 +40,50 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
		using VectorType = typename LocalFriction<dimension>::VectorType;
		using MatrixType = typename LocalFriction<dimension>::MatrixType;

		private:
		VectorType removeNormal(const VectorType& x) {
		VectorType res = x;
		res[0] = 0.0;
		return res;
		}

		public:
		template <typename... Args>
		WrappedScalarFriction(Args... args)
		: func_(args...) {}

		double operator()(VectorType const &x) const override {
		auto tangential_x = removeNormal(x);
		return func_.evaluate(tangential_x.two_norm());
		}

		void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }

		double regularity(VectorType const &x) const override {
		double const xnorm = x.two_norm();
		if (xnorm <= 0.0)
		auto tangential_x = removeNormal(x);
		double const xnorm = tangential_x.two_norm();
		if (xnorm < 0.0)
		return std::numeric_limits<double>::infinity();

		return func_.regularity(xnorm);
		}

		double coefficientOfFriction(VectorType const &x) const override {
		return func_.coefficientOfFriction(x.two_norm());
		auto tangential_x = removeNormal(x);
		return func_.coefficientOfFriction(tangential_x.two_norm());
		}

		// directional subdifferential: at u on the line u + t*v
		// u and v are assumed to be non-zero
		void directionalSubDiff(VectorType const &x, VectorType const &v,
		Dune::Solvers::Interval<double> &D) const override {
		double const xnorm = x.two_norm();

		auto tangential_x = removeNormal(x);
		double const xnorm = tangential_x.two_norm();
		if (xnorm <= 0.0)
		D[0] = D[1] = func_.differential(0.0) * v.two_norm();
		else
		D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
		D[0] = D[1] = func_.differential(xnorm) * (tangential_x * v) / xnorm;
		}

		/** Formula for the derivative:
		@@ -86,17 +103,26 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
		\f}
		*/
		void addHessian(VectorType const &x, MatrixType &A) const override {
		double const xnorm2 = x.two_norm2();
		//std::cout << A << std::endl;
		//std::cout << x << std::endl;

		/* double const xnorm2 = x.two_norm2();
		double const xnorm = std::sqrt(xnorm2);
		if (xnorm2 <= 0.0)
		return;

		//std::cout << xnorm << " " << xnorm2 << std::endl;

		double const H1 = func_.differential(xnorm);
		double const H2 = func_.second_deriv(xnorm);

		//std::cout << H1 << " " << H2 << std::endl;

		double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
		double const idweight = H1 / xnorm;

		//std::cout << tensorweight << " " << idweight << std::endl;

		for (size_t i = 0; i < dimension; ++i)
		for (size_t j = 0; j < i; ++j) {
		double const entry = tensorweight * x[i] * x[j];
		@@ -107,16 +133,47 @@ class WrappedScalarFriction : public LocalFriction<dimension> {
		for (size_t k = 0; k < dimension; ++k) {
		double const entry = tensorweight * x[k] * x[k];
		A[k][k] += entry + idweight;
		}*/
		auto tangential_x = removeNormal(x);
		const double xnorm = tangential_x.two_norm();

		if (std::isinf(func_.regularity(xnorm)) or xnorm<=0.0)
		return;

		VectorType y = tangential_x;
		y /= xnorm;

		double const H1 = func_.differential(xnorm);
		double const H2 = func_.second_deriv(xnorm);

		double const tensorweight = (H2 - H1 / xnorm);
		double const idweight = H1 / xnorm;

		//std::cout << tensorweight << " " << idweight << std::endl;

		for (size_t i = 0; i < dimension; ++i)
		for (size_t j = 0; j < i; ++j) {
		double const entry = tensorweight * y[i] * y[j];
		A[i][j] += entry;
		A[j][i] += entry;
		}

		for (size_t k = 0; k < dimension; ++k) {
		double const entry = tensorweight * y[k] * y[k];
		A[k][k] += entry + idweight;
		}

		//std::cout << A << std::endl;
		}

		void addGradient(VectorType const &x, VectorType &y) const override {
		double const xnorm = x.two_norm();
		auto tangential_x = removeNormal(x);
		double const xnorm = tangential_x.two_norm();
		if (std::isinf(func_.regularity(xnorm)))
		return;

		if (xnorm > 0.0)
		Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, x);
		Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, tangential_x);
		}

		void directionalDomain(VectorType const &, VectorType const &,

dune/tectonic/data-structures/globalfrictioncontainer.cc

0 → 100644

+114 −0

Original line number	Diff line number	Diff line
		#ifdef HAVE_CONFIG_H
		#include "config.h"
		#endif

		/*
		#include "globalfrictioncontainer.hh"

		template <class BaseGlobalFriction, size_t depth>
		auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i)
		-> IndexObject& {
		return globalFriction_[i];
		}

		template <class BaseGlobalFriction, size_t depth>
		auto GlobalFrictionContainer<BaseGlobalFriction, depth>::operator[](std::size_t i) const
		-> const IndexObject& {
		return globalFriction_[i];
		}

		template <class BaseGlobalFriction, size_t depth>
		auto GlobalFrictionContainer<BaseGlobalFriction, depth>::size() const
		-> size_t {
		return globalFriction_.size();
		}

		template <class BaseGlobalFriction, size_t depth>
		void GlobalFrictionContainer<BaseGlobalFriction, depth>::resize(std::list<size_t> list) {
		assert(list.size() <= depth);

		if (list.size() == 0) {
		globalFriction_.resize(0);
		} else {
		globalFriction_.resize(list.front());
		list.pop_front();

		for (size_t i=0; i<size(); i++) {
		globalFriction_[i].resize(list);
		}
		}
		}

		template <class BaseGlobalFriction, size_t depth>
		template <class VectorContainer>
		void GlobalFrictionContainer<BaseGlobalFriction, depth>::updateAlpha(const VectorContainer& newAlpha) {
		assert(newAlpha.size() == size());

		for (size_t i=0; i<size(); i++) {
		globalFriction_[i].updateAlpha(newAlpha[i]);
		}
		}

		template <class BaseGlobalFriction, size_t depth>
		auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction()
		-> GlobalFriction& {
		return globalFriction_;
		}

		template <class BaseGlobalFriction, size_t depth>
		auto GlobalFrictionContainer<BaseGlobalFriction, depth>::globalFriction() const
		-> const GlobalFriction& {
		return globalFriction_;
		}


		template <class BaseGlobalFriction>
		auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i)
		-> IndexObject& {
		return globalFriction_[i];
		}

		template <class BaseGlobalFriction>
		auto GlobalFrictionContainer<BaseGlobalFriction, 1>::operator[](std::size_t i) const
		-> const IndexObject& {
		return globalFriction_[i];
		}

		template <class BaseGlobalFriction>
		auto GlobalFrictionContainer<BaseGlobalFriction, 1>::size() const
		-> size_t {
		return globalFriction_.size();
		}

		template <class BaseGlobalFriction>
		void GlobalFrictionContainer<BaseGlobalFriction, 1>::resize(std::list<size_t> newSize) {
		if (newSize.size() > 0) {
		globalFriction_.resize(newSize.front(), nullptr);
		} else {
		globalFriction_.resize(0);
		}
		}

		template <class BaseGlobalFriction>
		template <class Vector>
		void GlobalFrictionContainer<BaseGlobalFriction, 1>::updateAlpha(const Vector& newAlpha) {
		for (size_t i=0; i<size(); i++) {
		globalFriction_[i]->updateAlpha(newAlpha);
		}
		}

		template <class BaseGlobalFriction>
		auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction()
		-> GlobalFriction& {
		return globalFriction_;
		}

		template <class BaseGlobalFriction>
		auto GlobalFrictionContainer<BaseGlobalFriction, 1>::globalFriction() const
		-> const GlobalFriction& {
		return globalFriction_;
		}


		#include "globalfrictioncontainer_tmpl.cc"
		*/

dune/tectonic/data-structures/globalfrictioncontainer.hh

0 → 100644

+114 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/matrices.hh

0 → 100644

+37 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/CMakeLists.txt

0 → 100644

+12 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/contactnetwork.cc

0 → 100644

+442 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/contactnetwork.hh

0 → 100644

+147 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/contactnetwork_tmpl.cc

0 → 100644

+10 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/levelcontactnetwork.cc

0 → 100644

+247 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/levelcontactnetwork.hh

0 → 100644

+111 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/network/levelcontactnetwork_tmpl.cc

0 → 100644

+14 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/data-structures/program_state.hh

0 → 100644

+245 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/explicitgrid.hh

0 → 100644

+18 −0

File added.

Preview size limit exceeded, changes collapsed.

src/explicitvectors.hh→dune/tectonic/explicitvectors.hh

+20 −0

File changed and moved.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/CMakeLists.txt

0 → 100644

+25 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/cantorfactory.cc

0 → 100644

+404 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/cantorfactory.hh

0 → 100644

+83 −0

File added.

Preview size limit exceeded, changes collapsed.

src/assemblers_tmpl.cc→dune/tectonic/factories/cantorfactory_tmpl.cc

+7 −0

File changed and moved.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/contactnetworkfactory.hh

0 → 100644

+59 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/levelcontactnetworkfactory.hh

0 → 100644

+54 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/stackedblocksfactory.cc

0 → 100644

+241 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/stackedblocksfactory.hh

0 → 100644

+87 −0

File added.

Preview size limit exceeded, changes collapsed.

src/hdf5/surface-writer_tmpl.cc→dune/tectonic/factories/stackedblocksfactory_tmpl.cc

+8 −0

File changed and moved.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/strikeslipfactory.cc

0 → 100644

+276 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/strikeslipfactory.hh

0 → 100644

+87 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/strikeslipfactory_tmpl.cc

0 → 100644

+8 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/threeblocksfactory.cc

0 → 100644

+204 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/threeblocksfactory.hh

0 → 100644

+87 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/threeblocksfactory_tmpl.cc

0 → 100644

+8 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/twoblocksfactory.cc

0 → 100644

+201 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/twoblocksfactory.hh

0 → 100644

+87 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/factories/twoblocksfactory_tmpl.cc

0 → 100644

+8 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/globalratestatefriction.hh

deleted100644 → 0

+0 −71

File deleted.

Preview size limit exceeded, changes collapsed.

src/gridselector.hh→dune/tectonic/gridselector.hh

+0 −0

File moved.

dune/tectonic/io/CMakeLists.txt

0 → 100644

+19 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5-bodywriter.hh

0 → 100644

+74 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5-writer.hh

0 → 100644

+111 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/CMakeLists.txt

0 → 100644

+25 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/frictionalboundary-writer.hh

0 → 100644

+86 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/iteration-writer.hh

0 → 100644

+42 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/restart-io.hh

0 → 100644

+62 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/restartbody-io.hh

0 → 100644

+52 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/restrict.hh

0 → 100644

+23 −0

File added.

Preview size limit exceeded, changes collapsed.

src/hdf5/surface-writer.cc→dune/tectonic/io/hdf5/surface-writer.hh

+53 −0

File changed and moved.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/hdf5/time-writer.hh

0 → 100644

+27 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/io-handler.hh

0 → 100644

+182 −0

File added.

Preview size limit exceeded, changes collapsed.

src/uniform-grid-writer.cc→dune/tectonic/io/uniform-grid-writer.cc

+7 −4

File changed and moved.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/vtk-bodywriter.hh

0 → 100644

+58 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/io/vtk.hh

0 → 100644

+49 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/minimisation.hh

deleted100644 → 0

+0 −41

File deleted.

Preview size limit exceeded, changes collapsed.

dune/tectonic/multi-threading/CMakeLists.txt

0 → 100644

+11 −0

File added.

Preview size limit exceeded, changes collapsed.

dune/tectonic/multi-threading/stoppable.hh

0 → 100644

+0 −0

Original line number	Diff line number	Diff line

Source

Target

Files

Some changes are not shown.