Make session less FE-centric

and more focused on interpolation (people might not be familar with FE).
issue#2
Tim Skov Jacobsen 2019-11-08 14:33:44 +01:00
parent e0bf1025f1
commit c5de8c7c8a
4 changed files with 500 additions and 242 deletions

View File

@ -287,39 +287,82 @@
"source": [
"# 8. Interpolation\n",
"\n",
"If you haven't already installed the packages `pandas`, `xlrd` and `scipy`, please do so.\n",
"*If you haven't already installed the packages `numpy`, `pandas`, `xlrd` and `scipy`, please do so.*\n",
"\n",
"This session has no new material per se. It's supposed to combine elements from the previous sessions into a larger exercise. \n",
"\n",
"The exercise is about 3D interpolation. A set of known points $(x_{known}, x_{known}, z_{known})$ have prescribed values and are used as basis for interpolating $z$-values for a large set of points where only $(x, y)$ are known.\n",
"\n",
"For performing the actual interpolation, we call a function from a third party library called `scipy`. It's built on top of of `numpy` and holds many operations used in scientific analysis. \n",
"\n",
"The code originates from a COWI project where the points $(x, y)$ represent node coordinates from a base slab in a Finite Element model, while $z$-coordinates denote settlement values. The known points $(x_{known}, y_{known}, z_{known})$ stem from a detailed geotechnical analysis which could only be performed in a certain amount of points. The settlement values in the remaining points $(x, y)$ were therefore put into the FE-model as imposed displacements by a procedure similar to this. \n",
"\n",
"This exercise is mostly about reading and understanding code written by others. Which is often just as important as being able to write it oneself.\n",
"\n",
"# Exercise 1.1\n",
"Read through the code given in the script below and try to understand what it does and how it does it.\n",
"\n",
"Copy the script to your editor and run it. \n",
"\n",
"Add print statements if you are unsure about how a certain variable looks at any point. Remember you can print the first five rows of a DataFrame with `df.head()`. \n",
"Add print statements if you are unsure about how a certain variable looks at any point throughout the code. Remember you can print the first five rows of a DataFrame with `df.head()`. \n",
"\n",
"The script reads two Excel files. One limitation of this is that they cannot be read while they are open. If you want to inspect these files while running the script, create a copy. \n",
"The script reads two Excel files, i.e. one containing known points and one containing points to be interpolated. \n",
"One limitation of Excel files is that they cannot be read while they are open. If you want to inspect these files while running the script, create a copy. \n",
"\n",
"# Exercise 1.2\n",
"The variable `settlements_interpolated` is a numpy array (you can see this by `type(settlements_interpolated)` => `<class 'numpy.ndarray'>`).\n",
"The bulk of the computational work in the script is done by the line:\n",
"\n",
"The last part of the code creates a figure object and an axis object which enables 3D plots.\n",
"---\n",
"```python\n",
"settlements_interpolated = griddata(xy_known, settlements_known, (x_nodes, y_nodes), method='cubic')\n",
"```\n",
"---\n",
"This is the `scipy` function that performs the interpolation. Try to read the documentation for it [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html)\n",
"\n",
"Continue the plotting code to add:\n",
"The returned value that we save in the variable `settlements_interpolated` is a numpy array. You can see this by `print(type(settlements_interpolated))` which returns: `<class 'numpy.ndarray'>`. \n",
"\n",
"* 3D scatter points of the known points $(x_{known}, y_{known}, w_{known})$\n",
"The last part of the given code creates a figure object and an axis object which enables 3D plots.\n",
"\n",
"* 3D scatter points of the interpolated settlements in every base slab node $(x_{nodes}, y_{nodes}, w_{nodes})$\n",
"**Continue the plotting code to add:**\n",
"\n",
"In the above $w$ denotes the settlements.\n",
"* 3D scatter plot of the known points $(x_{known}, y_{known}, z_{known})$\n",
"\n",
"* 3D scatter plot of the interpolated points $(x, y, z)$\n",
"\n",
"The plots should all be in the same figure and axis. Adding a scatter plot to an axis `ax` can be done as `ax.scatter(...)`.\n",
"\n",
"# Exercise 1.3\n",
"All the settlement values in the base slab are in this example to be transferred to and applied in a Sofistik calculation.\n",
"\n",
"To this end, we can create a text file with the exact syntax that the Sofistik input language accepts (this language is called CADINP and is somewhat similar to IBDAS's input language). The file has to be a `.dat`-file.\n",
"As mentioned, this was used in a project for interpolating settlement values to be applied in an FE-model. The FE-software (Sofistik) has a certain input language which the interpolated values needed to be blended into. \n",
"\n",
"To write data to a file we can use something called **context managers**. Basically, it allows us to open a file and write to it. See code snippet below:\n",
"In this exercise we will construct the input `.dat`.file that Sofistik can read. \n",
"\n",
">**Note:** This procedure could be used to produce many file types. It's a good way to programmatically create input files to software. This is just one specific example of a use case.\n",
"\n",
"The Sofistik input file we want to create has this syntax:\n",
"\n",
"---\n",
"\n",
"<pre><code>\n",
"+PROG SOFILOAD \n",
"\n",
"LC 25 type 'SL' fact 1.0 facd 0.0 titl 'LT settlement all nodes'\n",
"\n",
" POIN NODE <font color='#1E90FF'>insert first node number</font> WIDE 0 TYPE WZZ <font color='#1E90FF'>insert first interpolated z-value</font>\n",
" ...\n",
" <font color='#1E90FF'><i>one line per pair of node number/z-value</i></font> \n",
" ...\n",
" POIN NODE <font color='#1E90FF'>insert last node number</font> WIDE 0 TYPE WZZ <font color='#1E90FF'>insert last interpolated z-value</font>\n",
" \n",
"END\n",
"</code></pre>\n",
"\n",
"---\n",
"\n",
"The indented block should print all the node/settlement pairs. The three non-indented lines should only appear once. The output file should look like the file `interpolation_output_example.dat` in the folder. Newlines are made by `\\n`.\n",
"\n",
"\n",
"## How to write to files\n",
"To write data to a file we can use something called a **context manager**. Basically, it allows us to open a file and write to it. See code snippet below:\n",
"\n",
"---\n",
"~~~python\n",
@ -331,30 +374,14 @@
"~~~\n",
"---\n",
"\n",
"By using the concept the file is automatically closed after our indented block is terminated. It also creates the file in case it doesn't already exist. \n",
"\n",
"The Sofistik input file we want has the format:\n",
"\n",
"~~~html\n",
"+PROG SOFILOAD \n",
"\n",
"LC 25 type 'P' fact 1.0 facd 0.0 titl 'LT settlement all nodes'\n",
"\n",
" POIN NODE 'insert first node number' WIDE 0 TYPE WZZ 'insert first settlement value'\n",
" ... 'one line per base slab node' ... \n",
" POIN NODE 'insert last node number' WIDE 0 TYPE WZZ 'insert last settlement value'\n",
" \n",
"END\n",
"~~~\n",
"\n",
"The indented block should print all the node/settlement pairs. The three non-indented lines should only appear once. The output file should look like the file `settlement_field_generated_code_example.dat` in the folder.\n"
"By using the concept the file is automatically closed after our indented block is terminated. It also creates the file in case it doesn't already exist. \n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The script\n",
"# The script\n",
"\n",
"~~~python\n",
"\n",
@ -365,24 +392,24 @@
"from mpl_toolkits.mplot3d import Axes3D\n",
"\n",
"\n",
"# Set name of Excel file to read\n",
"file_known = 'known_sections_plaxis.xlsx'\n",
"# Set name of Excel file to read containing known points\n",
"file_known = 'known_points.xlsx'\n",
"\n",
"# Set name of sheet to read from Excel file\n",
"sheet_known = 'known_sections_plaxis'\n",
"sheet_known = 'Sheet1'\n",
"\n",
"# Read data from Excel sheet into a dataframe\n",
"df = pd.read_excel(file_known, sheet_name=sheet_known, skiprows=7)\n",
"\n",
"# Extract columns whose names starts with a 'Y' into new dataframe of Y-coordinates\n",
"# Extract column names starting with 'Y' into new dataframe of known Y-coords\n",
"df_y = df[df.columns[df.columns.str.startswith('Y')]]\n",
"\n",
"# Extract columns whose names starts with 'Z' into new dataframe of settlements\n",
"df_settlements_known = df[df.columns[df.columns.str.startswith('Z')]]\n",
"# Extract column names starting with 'Z' into new dataframe of known Z-coords\n",
"df_z_known = df[df.columns[df.columns.str.startswith('Z')]]\n",
"\n",
"# Flatten dataframe values into 1D array\n",
"# Flatten dataframe values into 1D array (matri format -> vector format)\n",
"y_known = df_y.values.flatten()\n",
"settlements_known = df_settlements_known.values.flatten()\n",
"z_known = df_z_known.values.flatten()\n",
"\n",
"# Extract known x-values\n",
"x_known = df['X']\n",
@ -391,28 +418,28 @@
"# This will create matching(x, y)-points between arrays x and y\n",
"x_known = np.repeat(x_known, len(df_y.columns))\n",
"\n",
"# Set names and read Excel file with base slab nodes\n",
"file_nodes = 'base_slab_nodes.xlsx'\n",
"sheet_nodes = 'XLSX-Export'\n",
"df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)\n",
"\n",
"# Extract x- and y-coordinates of nodes\n",
"x_nodes = df_nodes['X [m]']\n",
"y_nodes = df_nodes['Y [m]']\n",
"\n",
"# Extract node numbers\n",
"node_no = df_nodes['NR']\n",
"\n",
"# Mirror known y-values and add corresponding x-values and settlements\n",
"# Mirror known y-values and add corresponding x- and y-values\n",
"x_known = np.append(x_known, x_known)\n",
"y_known = np.append(y_known, -y_known)\n",
"settlements_known = np.append(settlements_known, settlements_known)\n",
"z_known = np.append(z_known, z_known)\n",
"\n",
"# Arrange known (x, y) points to fit input for interpolation\n",
"xy_known = np.array(list(zip(x_known, y_known)))\n",
"\n",
"# Set names and read Excel file with nodes to be interpolated\n",
"file_nodes = 'points_to_be_interpolated.xlsx'\n",
"sheet_nodes = 'XLSX-Export'\n",
"df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)\n",
"\n",
"# Extract x- and y-coordinates of nodes to be interpolated\n",
"x_nodes = df_nodes['X [m]']\n",
"y_nodes = df_nodes['Y [m]']\n",
"\n",
"# Extract node numbers for points to be interpolated\n",
"node_no = df_nodes['NR']\n",
"\n",
"# Perform interpolation calculation\n",
"settlements_interpolated = griddata(xy_known, settlements_known, (x_nodes, y_nodes), method='cubic')\n",
"points_interpolated = griddata(xy_known, z_known, (x_nodes, y_nodes), method='cubic')\n",
"\n",
"\n",
"####################\n",
@ -424,10 +451,10 @@
"# Create axis object for 3D plot\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"\n",
"# Plot known settlement points as 3D scatter plot (ax.scatter(...))\n",
"# Plot known points as 3D scatter plot (ax.scatter(...))\n",
" # <Put plotting code here!>\n",
"\n",
"# Plot interpolated field as 3D scatter plot\n",
"# Plot interpolated points as 3D scatter plot\n",
" # <Put plotting code here!>\n",
"\n",
"# Show figure\n",
@ -437,10 +464,9 @@
"####################\n",
"### Exercise 1.3 ###\n",
"####################\n",
"# Write Sofistik input code to .dat-file for applying settlements as imposed displacements\n",
" # <Put plotting code here!>\n",
" \n",
" \n",
"# Write Sofistik input code to .dat-file for applying the interpolated z-values as \n",
"# imposed displacement load in all points (x, y)\n",
" # <Put code that creates and writes to a .dat file here!> \n",
"\n",
"~~~"
]

View File

@ -1,10 +1,291 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<style>\n",
"/* div#notebook {\n",
" font-size: 13pt;\n",
" line-height: 120%;\n",
" color: #303030;\n",
" -webkit-font-smoothing: antialiased !important;\n",
" padding-top: 25px !important;\n",
"} */\n",
"\n",
"\n",
"body,\n",
"div.body {\n",
" font-family: Roboto;\n",
" font-size: 16pt;\n",
" color: #303030;\n",
" /* background-color: #d1b5b5; */\n",
" /* background: #8f4f4f; */\n",
" margin-right: 10px;\n",
" -webkit-font-smoothing: antialiased !important;\n",
"}\n",
"\n",
"/* Code inside HTML/Markdown */\n",
"div.rendered_html code {\n",
" border-radius: 5px;\n",
"}\n",
"\n",
"/* Output area from code cells */\n",
"div.output_area pre {\n",
" font-size: 11pt !important;\n",
" color: #303030;\n",
"}\n",
"\n",
"\n",
"\n",
"/* HEADING 1 styles */\n",
"h1 {\n",
" font-family: 'Roboto', 'Signika', sans-serif;\n",
" font-weight: ligher;\n",
" font-style: normal;\n",
" font-size: 20pt;\n",
" width: 100%;\n",
" text-align: left;\n",
" color: #EE7631;\n",
" border-bottom: 1px solid rgba(238, 118, 49, 0.575);\n",
" padding-bottom: 6px;\n",
" /* font-variant: small-caps; */\n",
" }\n",
" \n",
" /* table {\n",
" margin: 5px;\n",
" width: 290px;\n",
" }\n",
" \n",
" th {\n",
" padding: 3px;\n",
" }\n",
" \n",
" td {\n",
" padding-left: 8px;\n",
" padding-right: 8px;\n",
" border: 1px solid #990000;\n",
" background-color: #ffffcc;\n",
" }\n",
"\n",
" #trHeader {\n",
" text-decoration: underline;\n",
" color: #990000;\n",
" }\n",
" \n",
" .centerCell {\n",
" text-align: center;\n",
" } */\n",
"\n",
"/* HEADING 2 styles */\n",
"h2 {\n",
" font-family: \"Roboto\";\n",
" text-align: left;\n",
" font-size: 12pt;\n",
" color: #EE7631;\n",
" margin-bottom: 5px;\n",
" /* border-bottom: 1px solid lightgray; */\n",
" border-bottom: 0.8px solid rgba(238, 118, 49, 0.575);\n",
" padding-bottom: 6px\n",
" }\n",
" \n",
" /* table {\n",
" margin: 5px;\n",
" width: 290px;\n",
" }\n",
" \n",
" th {\n",
" padding: 3px;\n",
" }\n",
" \n",
" td {\n",
" padding-left: 8px;\n",
" padding-right: 8px;\n",
" border: 1px solid #990000;\n",
" background-color: #ffffcc;\n",
" }\n",
"\n",
" #trHeader {\n",
" text-decoration: underline;\n",
" color: #990000;\n",
" }\n",
" \n",
" .centerCell {\n",
" text-align: center;\n",
" } */\n",
"\n",
"h3 {\n",
" font-family: 'Roboto';\n",
" /* text-align: left; */\n",
" /* font-size: 12pt; */\n",
" /* color: #EE7631; */\n",
" /* margin-bottom: 5px; */\n",
" border-bottom: 0.5px solid #ededed;\n",
" padding-bottom: 6px\n",
" }\n",
"\n",
"\n",
"\n",
"p {\n",
" font-family: Roboto;\n",
" font-size: 16px;\n",
"}\n",
"\n",
"/* Lists with dots */\n",
"ul {\n",
" font-size: 16px;\n",
" line-height: 150%;\n",
"}\n",
"\n",
"/* Lists with numbers */\n",
"ol {\n",
" font-size: 16px;\n",
" line-height: 150%;\n",
"}\n",
"\n",
"/* Horizontal rules */\n",
"hr { \n",
" margin-top: 3px; \n",
" margin-bottom: 3px \n",
"}\n",
"\n",
"/* Links */\n",
"a {\n",
" color: #EE7631;\n",
"}\n",
"\n",
"\n",
"/* Code cells in the notebook - NOTE: color is font-color */\n",
".cm-s-ipython.CodeMirror {\n",
" font-family: monospace, monospace;\n",
" font-size: 11pt;\n",
" background: #ededed;\n",
" color: #303030; \n",
" border-radius: 2px;\n",
" /* margin-right: 10px; */\n",
" font-style: normal;\n",
" font-weight: normal;\n",
"}\n",
"\n",
"\n",
"/* Background of code cells */\n",
".cm-s-ipython.CodeMirror {\n",
" font-family: monospace, monospace;\n",
" font-size: 11pt;\n",
" background: rgba(211, 211, 211, 0.123);\n",
" color: #303030;\n",
" border-radius: 2px;\n",
" font-style: normal;\n",
" font-weight: normal;\n",
"}\n",
"\n",
"\n",
"/* .CodeMirror-gutters {\n",
" border: none;\n",
" border-right: 1px solid #e0e1e3 !important;\n",
" background-color: #e0e1e3 !important;\n",
" background: #e0e1e3 !important;\n",
" border-radius: 0px;\n",
" white-space: nowrap;\n",
"} */\n",
"\n",
"\n",
"\n",
"/* Code syntax highlithig theme */\n",
".cm-s-ipython .CodeMirror-cursor {\n",
" border-left: 2px solid #ff711a !important;\n",
"}\n",
".cm-s-ipython span.cm-comment {\n",
" color: #8d8d8d;\n",
" font-style: italic; \n",
"}\n",
".cm-s-ipython span.cm-atom {\n",
" color: #055be0;\n",
"}\n",
".cm-s-ipython span.cm-number {\n",
" color: #ff8132;\n",
"}\n",
".cm-s-ipython span.cm-property {\n",
" color: #303030;\n",
"}\n",
".cm-s-ipython span.cm-attribute {\n",
" color: #303030;\n",
"}\n",
".cm-s-ipython span.cm-keyword {\n",
" color: #a045ddf3;\n",
" font-weight: normal;\n",
"}\n",
".cm-s-ipython span.cm-string {\n",
" color: #009e07;\n",
"}\n",
".cm-s-ipython span.cm-meta {\n",
" color: #aa22ff;\n",
"}\n",
".cm-s-ipython span.cm-operator {\n",
" color: #055be0;\n",
"}\n",
".cm-s-ipython span.cm-builtin {\n",
" color: #3f2ce7;\n",
"}\n",
".cm-s-ipython span.cm-variable {\n",
" color: #303030;\n",
"}\n",
".cm-s-ipython span.cm-variable-2 {\n",
" color: #de143d;\n",
"}\n",
".cm-s-ipython span.cm-variable-3 {\n",
" color: #aa22ff;\n",
"}\n",
".cm-s-ipython span.cm-def {\n",
" color: #e22978;\n",
" font-weight: normal;\n",
"}\n",
".cm-s-ipython span.cm-error {\n",
" background: rgba(191,11,55,.70);\n",
"}\n",
".cm-s-ipython span.cm-tag {\n",
" color: #e22978;\n",
"}\n",
".cm-s-ipython span.cm-link {\n",
" color: #ef5c00;\n",
"}\n",
".cm-s-ipython span.cm-storage {\n",
" color: #055be0;\n",
"}\n",
".cm-s-ipython span.cm-entity {\n",
" color: #e22978;\n",
"}\n",
".cm-s-ipython span.cm-quote {\n",
" color: #009e07;\n",
"}\n",
"\n",
"\n",
"\n",
"</style>"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import HTML\n",
"HTML('<style>{}</style>'.format(open('../css/cowi.css').read()))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Session 8 - Exercise solution\n",
"# 8. Exercise solution\n",
"\n",
"The full script is provided below."
]
@ -21,24 +302,24 @@
"from mpl_toolkits.mplot3d import Axes3D\n",
"\n",
"\n",
"# Set name of Excel file to read\n",
"file_known = 'known_sections_plaxis.xlsx'\n",
"# Set name of Excel file to read containing known points\n",
"file_known = 'known_points.xlsx'\n",
"\n",
"# Set name of sheet to read from Excel file\n",
"sheet_known = 'known_sections_plaxis'\n",
"sheet_known = 'Sheet1'\n",
"\n",
"# Read data from Excel sheet into a dataframe\n",
"df = pd.read_excel(file_known, sheet_name=sheet_known, skiprows=7)\n",
"\n",
"# Extract columns whose names starts with a 'Y' into new dataframe of Y-coordinates\n",
"# Extract column names starting with 'Y' into new dataframe of known Y-coords\n",
"df_y = df[df.columns[df.columns.str.startswith('Y')]]\n",
"\n",
"# Extract columns whose names starts with 'Z' into new dataframe of settlements\n",
"df_settlements_known = df[df.columns[df.columns.str.startswith('Z')]]\n",
"# Extract column names starting with 'Z' into new dataframe of known Z-coords\n",
"df_z_known = df[df.columns[df.columns.str.startswith('Z')]]\n",
"\n",
"# Flatten dataframe values into 1D array\n",
"# Flatten dataframe values into 1D array (matri format -> vector format)\n",
"y_known = df_y.values.flatten()\n",
"settlements_known = df_settlements_known.values.flatten()\n",
"z_known = df_z_known.values.flatten()\n",
"\n",
"# Extract known x-values\n",
"x_known = df['X']\n",
@ -47,28 +328,28 @@
"# This will create matching(x, y)-points between arrays x and y\n",
"x_known = np.repeat(x_known, len(df_y.columns))\n",
"\n",
"# Set names and read Excel file with base slab nodes\n",
"file_nodes = 'base_slab_nodes.xlsx'\n",
"sheet_nodes = 'XLSX-Export'\n",
"df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)\n",
"\n",
"# Extract x- and y-coordinates of nodes\n",
"x_nodes = df_nodes['X [m]']\n",
"y_nodes = df_nodes['Y [m]']\n",
"\n",
"# Extract node numbers\n",
"node_no = df_nodes['NR']\n",
"\n",
"# Mirror known y-values and add corresponding x-values and settlements\n",
"# Mirror known y-values and add corresponding x- and y-values\n",
"x_known = np.append(x_known, x_known)\n",
"y_known = np.append(y_known, -y_known)\n",
"settlements_known = np.append(settlements_known, settlements_known)\n",
"z_known = np.append(z_known, z_known)\n",
"\n",
"# Arrange known (x, y) points to fit input for interpolation\n",
"xy_known = np.array(list(zip(x_known, y_known)))\n",
"\n",
"# Set names and read Excel file with nodes to be interpolated\n",
"file_nodes = 'points_to_be_interpolated.xlsx'\n",
"sheet_nodes = 'XLSX-Export'\n",
"df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)\n",
"\n",
"# Extract x- and y-coordinates of nodes to be interpolated\n",
"x_nodes = df_nodes['X [m]']\n",
"y_nodes = df_nodes['Y [m]']\n",
"\n",
"# Extract node numbers for points to be interpolated\n",
"node_no = df_nodes['NR']\n",
"\n",
"# Perform interpolation calculation\n",
"settlements_interpolated = griddata(xy_known, settlements_known, (x_nodes, y_nodes), method='cubic')\n",
"z_interpolated = griddata(xy_known, z_known, (x_nodes, y_nodes), method='cubic')\n",
"\n",
"\n",
"####################\n",
@ -77,14 +358,15 @@
"# Create figure object\n",
"fig = plt.figure()\n",
"\n",
"# Create axis object for 3D plot and put it on the figure\n",
"# Create axis object for 3D plot\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"\n",
"# Plot known settlement points as 3D scatter plot (ax.scatter(...))\n",
"ax.scatter(x_known, y_known, settlement_known, '-.', color='limegreen')\n",
"# Plot known points as 3D scatter plot (ax.scatter(...))\n",
"ax.scatter(x_known, y_known, z_known, '-.', color='limegreen')\n",
"\n",
"# Plot interpolated field as 3D scatter plot\n",
"ax.scatter(x_nodes, y_nodes, settlement_interpolated, '.', color='cornflowerblue', s=0.1)\n",
"# Plot interpolated points as 3D scatter plot\n",
"ax.scatter(x_nodes, y_nodes, z_interpolated,\n",
" '.', color='cornflowerblue', s=0.1)\n",
"\n",
"# Show figure\n",
"plt.show()\n",
@ -93,18 +375,19 @@
"####################\n",
"### Exercise 1.3 ###\n",
"####################\n",
"# Write Sofistik input code to .dat-file for applying settlements as imposed displacements\n",
"with open(f'Generated_sofistik_code.dat', 'w') as file:\n",
" \n",
" # Write the 'static' text to file \n",
"# Write Sofistik input code to .dat-file for applying the interpolated z-values as\n",
"# imposed displacement load (settlement) in all points (x, y)\n",
"with open(f'generated_file.dat', 'w') as file:\n",
"\n",
" # Write the 'static' text to file\n",
" file.write('''+PROG SOFILOAD \n",
"\n",
"LC 25 type 'P' fact 1.0 facd 0.0 titl 'LT settlement all nodes' \\n''')\n",
" \n",
"\n",
" # Write the 'variable' text to file with node number/settlement pairs\n",
" for node, settlement in zip(node_no, settlement_interpolated):\n",
" for node, settlement in zip(node_no, z_interpolated):\n",
" file.write(f' POIN NODE {node} WIDE 0 TYPE WZZ {settlement} \\n')\n",
" \n",
"\n",
" # Write 'static' END statement to file\n",
" file.write('END')\n",
"\n",
@ -115,7 +398,7 @@
"metadata": {
"hide_input": false,
"kernelspec": {
"display_name": "Python [default]",
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
@ -129,7 +412,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
"version": "3.7.5"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,

View File

@ -5,7 +5,7 @@ import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Set name of Excel file to read
# Set name of Excel file to read containing known points
file_known = 'known_points.xlsx'
# Set name of sheet to read from Excel file
@ -14,15 +14,15 @@ sheet_known = 'Sheet1'
# Read data from Excel sheet into a dataframe
df = pd.read_excel(file_known, sheet_name=sheet_known, skiprows=7)
# Extract columns whose names starts with a 'Y' into new dataframe of Y-coordinates
# Extract column names starting with 'Y' into new dataframe of known Y-coords
df_y = df[df.columns[df.columns.str.startswith('Y')]]
# Extract columns whose names starts with 'Z' into new dataframe of settlements
df_settlements_known = df[df.columns[df.columns.str.startswith('Z')]]
# Extract column names starting with 'Z' into new dataframe of known Z-coords
df_z_known = df[df.columns[df.columns.str.startswith('Z')]]
# Flatten dataframe values into 1D array
# Flatten dataframe values into 1D array (matri format -> vector format)
y_known = df_y.values.flatten()
settlements_known = df_settlements_known.values.flatten()
z_known = df_z_known.values.flatten()
# Extract known x-values
x_known = df['X']
@ -31,28 +31,28 @@ x_known = df['X']
# This will create matching(x, y)-points between arrays x and y
x_known = np.repeat(x_known, len(df_y.columns))
# Set names and read Excel file with base slab nodes
file_nodes = 'points_to_be_interpolated.xlsx'
sheet_nodes = 'XLSX-Export'
df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)
# Extract x- and y-coordinates of nodes
x_nodes = df_nodes['X [m]']
y_nodes = df_nodes['Y [m]']
# Extract node numbers
node_no = df_nodes['NR']
# Mirror known y-values and add corresponding x-values and settlements
# Mirror known y-values and add corresponding x- and y-values
x_known = np.append(x_known, x_known)
y_known = np.append(y_known, -y_known)
settlements_known = np.append(settlements_known, settlements_known)
z_known = np.append(z_known, z_known)
# Arrange known (x, y) points to fit input for interpolation
xy_known = np.array(list(zip(x_known, y_known)))
# Set names and read Excel file with nodes to be interpolated
file_nodes = 'points_to_be_interpolated.xlsx'
sheet_nodes = 'XLSX-Export'
df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)
# Extract x- and y-coordinates of nodes to be interpolated
x_nodes = df_nodes['X [m]']
y_nodes = df_nodes['Y [m]']
# Extract node numbers for points to be interpolated
node_no = df_nodes['NR']
# Perform interpolation calculation
settlements_interpolated = griddata(xy_known, settlements_known, (x_nodes, y_nodes), method='cubic')
points_interpolated = griddata(xy_known, z_known, (x_nodes, y_nodes), method='cubic')
####################
@ -64,10 +64,10 @@ fig = plt.figure()
# Create axis object for 3D plot
ax = fig.add_subplot(111, projection='3d')
# Plot known settlement points as 3D scatter plot (ax.scatter(...))
# Plot known points as 3D scatter plot (ax.scatter(...))
# <Put plotting code here!>
# Plot interpolated field as 3D scatter plot
# Plot interpolated points as 3D scatter plot
# <Put plotting code here!>
# Show figure
@ -77,6 +77,7 @@ ax = fig.add_subplot(111, projection='3d')
####################
### Exercise 1.3 ###
####################
# Write Sofistik input code to .dat-file for applying settlements as imposed displacements
# <Put plotting code here!>
# Write Sofistik input code to .dat-file for applying the interpolated z-values as
# imposed displacement load in all points (x, y)
# <Put code that creates and writes to a .dat file here!>

View File

@ -1,147 +1,95 @@
import numpy as np
import pandas as pd
import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
### SET VERSION NUMBER ###
version_number = 'v1a'
### SET WHAT TO PLOT ###
plot_model_points = False
plot_known_points = True
plot_interpolated_points = True
# Set name of Excel file to read containing known points
file_known = 'known_points.xlsx'
### READ EXCEL FILE WITH SECTIONS AND KNOWN SETTLEMENT DATA ###
file_name = f'known_sections_plaxis_{version_number}.xlsx'
sheet_name = 'known_sections_plaxis'
df_x = pd.read_excel(file_name, sheet_name=sheet_name, skiprows=7, usecols=[1])
df_y = pd.read_excel(file_name, sheet_name=sheet_name, skiprows=7, usecols=[2,3,4,5])
df_z25 = pd.read_excel(file_name, sheet_name=sheet_name, skiprows=7, usecols=[6,7,8,9])
df_z26 = pd.read_excel(file_name, sheet_name=sheet_name, skiprows=7, usecols=[10,11,12,13])
df_z27 = pd.read_excel(file_name, sheet_name=sheet_name, skiprows=7, usecols=[14,15,16,17])
# Set name of sheet to read from Excel file
sheet_known = 'Sheet1'
# Insert as many X-colunms as there are columns in the other dataframes
for i in range(2, len(df_y.columns)+1):
df_x[f'X{i}'] = df_x['X']
# Read data from Excel sheet into a dataframe
df = pd.read_excel(file_known, sheet_name=sheet_known, skiprows=7)
# Flatten dataframe of y-values into array
y_temp = df_y.values.flatten()
# Extract column names starting with 'Y' into new dataframe of known Y-coords
df_y = df[df.columns[df.columns.str.startswith('Y')]]
# Get indecies where values are not NaN
idx = np.where(~np.isnan(y_temp))
# Extract column names starting with 'Z' into new dataframe of known Z-coords
df_z_known = df[df.columns[df.columns.str.startswith('Z')]]
# Extract all non-NaN values to create array of all known points
y_known = y_temp[idx]
x_known = df_x.values.flatten()[idx]
z25_vals = df_z25.values.flatten()[idx]
z26_vals = df_z26.values.flatten()[idx]
z27_vals = df_z27.values.flatten()[idx]
# Flatten dataframe values into 1D array (matri format -> vector format)
y_known = df_y.values.flatten()
z_known = df_z_known.values.flatten()
# Extract known x-values
x_known = df['X']
# Mirror all data, since everything is assumed symmetric about CL (only vals defined on one side)
# Create X-array by repeating itself as many times as there are Y-columns
# This will create matching(x, y)-points between arrays x and y
x_known = np.repeat(x_known, len(df_y.columns))
# Mirror known y-values and add corresponding x- and y-values
x_known = np.append(x_known, x_known)
y_known = np.append(y_known, -y_known)
settlements_known = {'25': np.append(z25_vals, z25_vals),
'26': np.append(z26_vals, z26_vals),
'27': np.append(z27_vals, z27_vals)}
z_known = np.append(z_known, z_known)
for lc, settlement_known in settlements_known.items():
# Arrange known (x, y) points to fit input for interpolation
xy_known = np.array(list(zip(x_known, y_known)))
# FIXME A lot of stuff should go outside this loop!
# Set names and read Excel file with nodes to be interpolated
file_nodes = 'points_to_be_interpolated.xlsx'
sheet_nodes = 'XLSX-Export'
df_nodes = pd.read_excel(file_nodes, sheet_name=sheet_nodes)
### D-WALL NODES ###
# Read data file for structural points and their coordinates (To get D-wall bottom points)
dfp = pd.read_excel(f'structural_points_{version_number}.xlsx', sheet_name='XLSX-Export')
# Extract x- and y-coordinates of nodes to be interpolated
x_nodes = df_nodes['X [m]']
y_nodes = df_nodes['Y [m]']
# Remove leading or trailing white space from column names
dfp.columns = dfp.columns.str.strip()
# Extract node numbers for points to be interpolated
node_no = df_nodes['NR']
# Extract D-wall bottom corner points to new dataframe
df_dwalls_trumpet = dfp[(dfp['Text'] == 'Point') & (dfp['NR'] >= 808) & (dfp['NR'] <= 929)]
df_dwalls_station1 = dfp[(dfp['Text'] == 'Point') & (dfp['NR'] >= 1021) & (dfp['NR'] <= 1047)]
df_dwalls_station2 = dfp[(dfp['Text'] == 'Point') & (dfp['NR'] >= 1051) & (dfp['NR'] <= 1077)]
df_dwalls = pd.concat([df_dwalls_trumpet, df_dwalls_station1, df_dwalls_station2], axis=0)
# Gather corner points of D-walls (x, y, z) where settlement values will be interpolated
x_dwalls, y_dwalls, z_dwalls = df_dwalls['X [m]'], df_dwalls['Y [m]'], -df_dwalls['Z [m]']
node_no_dwalls = df_dwalls['NR']
### BASE SLAB NODES ###
# Read file with base slab node numbers and their coordinates into a dataframe
df_slabs = pd.read_excel(f'base_slab_nodes_{version_number}.xlsx', sheet_name='XLSX-Export')
# Remove leading or trailing white space from column names
df_slabs.columns = df_slabs.columns.str.strip()
x_slabs, y_slabs, z_slabs = df_slabs['X [m]'], df_slabs['Y [m]'], df_slabs['Z [m]']
node_no_slabs = df_slabs['NR']
### COMBINE BASE SLAB AND D-WALL DATA ###
x_nodes = np.append(x_dwalls, x_slabs)
y_nodes = np.append(y_dwalls, y_slabs)
z_nodes = np.append(z_dwalls, z_slabs)
node_no = np.append(node_no_dwalls, node_no_slabs)
### PERFORM INTERPOLATION ###
# x-y coordinates of points with known displacements
xy_known = np.array(list(zip(x_known, y_known))) # FIXME SHould be outside loop
# Calculate the interpolated z-values
settlement_interpolated = griddata(xy_known, settlement_known, (x_nodes, y_nodes), method='linear')
# Check if interpolated settlements have any nan values
print('-------------------------------------------')
print(f'LC{lc}:')
if np.isnan(settlement_interpolated).any():
print(' INFO:')
print(" Some interpolated settlement values are 'nan'.")
print(' This is probably beacause they fall out of the region defined by known points')
print(' (extrapolation not suported).')
nans = np.argwhere(np.isnan(settlement_interpolated))
print(f' Number of nan values are {len(nans)} out of {len(settlement_interpolated)} total values.')
print(f' X-values for points that failed to intepolate are:\n{x_nodes[nans.flatten()]}')
# TODO: Replace those nan-values with 0 and make sure they are colored red in the plot
else:
print(' All values interpolated succesfully!')
print('-------------------------------------------')
# Perform interpolation calculation
z_interpolated = griddata(xy_known, z_known, (x_nodes, y_nodes), method='cubic')
### WRITE INTERPOLATED FIELD TO .DAT FILE AS TEDDY CODE ###
# Write Teddy code for applying interpolated settlements to file
with open(f'teddy_code_settlement_field_LC{lc}_{version_number}.dat', 'w') as file:
# file.write('''+PROG SOFILOAD urs:47.8 $ Plaxis settlement
# HEAD Plaxis settlement-LT (interp. all nodes)
# UNIT TYPE 5
#
# LC 25 type 'P' fact 1.0 facd 0.0 titl 'LT settlement all nodes' \n''')
for node, settlement in zip(node_no, settlement_interpolated):
file.write(f' POIN NODE {node} WIDE 0 TYPE WZZ {settlement} \n')
# file.write('END')
####################
### Exercise 1.2 ###
####################
# Create figure object
fig = plt.figure()
# Create axis object for 3D plot
ax = fig.add_subplot(111, projection='3d')
# Plot known points as 3D scatter plot (ax.scatter(...))
ax.scatter(x_known, y_known, z_known, '-.', color='limegreen')
# Plot interpolated points as 3D scatter plot
ax.scatter(x_nodes, y_nodes, z_interpolated,
'.', color='cornflowerblue', s=0.1)
# Show figure
plt.show()
### PLOT KNOWN POINTS WITH INTERPOLATED FIELD ###
# Check of any plot options are set to True
if any([plot_model_points, plot_known_points, plot_interpolated_points]):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
####################
### Exercise 1.3 ###
####################
# Write Sofistik input code to .dat-file for applying the interpolated z-values as
# imposed displacement load (settlement) in all points (x, y)
with open(f'generated_file.dat', 'w') as file:
if plot_model_points:
# Plot structural points/nodes in the model
ax.scatter(x_nodes, y_nodes, z_nodes, '.', s=0.1) # Base slabs
ax.scatter(x_dwalls, y_dwalls, z_dwalls, '.', s=0.1, color='magenta') # D-walls
# Write the 'static' text to file
file.write('''+PROG SOFILOAD
if plot_known_points:
# Plot known settlement points
ax.scatter(x_known, y_known, settlement_known, '-.', color='limegreen')
LC 25 type 'SL' fact 1.0 facd 0.0 titl 'LT settlement all nodes' \n''')
if plot_interpolated_points:
# Plot interpolated field
ax.scatter(x_nodes, y_nodes, settlement_interpolated, '.', color='cornflowerblue', s=0.1)
# Write the 'variable' text to file with node number/settlement pairs
for node, settlement in zip(node_no, z_interpolated):
file.write(f' POIN NODE {node} WIDE 0 TYPE WZZ {settlement} \n')
# Set limits
# ax.set_xlim(6800, 7350)
# ax.set_zlim(-22, -15)
# ax.set_ylim(-100, 100)
plt.show()
# Write 'static' END statement to file
file.write('END')