From 309fc642c9299eca4e1665c069367c5cf5d18fef Mon Sep 17 00:00:00 2001 From: Krishna Kumar Date: Tue, 1 Oct 2024 12:50:56 -0500 Subject: [PATCH] OpenSees example --- examples/opensees/DS_GenFunctionsV3.py | 93 + examples/opensees/N10_T3.tcl | 515 ++ examples/opensees/ShortReport.tex | 98 + examples/opensees/interactiveplot.py | 73 + examples/opensees/opensees-dapi-v3.ipynb | 577 ++ examples/opensees/plotAcc.py | 64 + examples/opensees/plotPorepressure.py | 35 + examples/opensees/plotProfile.py | 87 + examples/opensees/plotStressStrain.py | 26 + examples/opensees/respSpectra.py | 70 + examples/opensees/schematic.png | Bin 0 -> 7794 bytes examples/opensees/short.tex | 158 + examples/opensees/velocity.input | 7998 ++++++++++++++++++++++ 13 files changed, 9794 insertions(+) create mode 100644 examples/opensees/DS_GenFunctionsV3.py create mode 100644 examples/opensees/N10_T3.tcl create mode 100644 examples/opensees/ShortReport.tex create mode 100644 examples/opensees/interactiveplot.py create mode 100644 examples/opensees/opensees-dapi-v3.ipynb create mode 100644 examples/opensees/plotAcc.py create mode 100644 examples/opensees/plotPorepressure.py create mode 100644 examples/opensees/plotProfile.py create mode 100644 examples/opensees/plotStressStrain.py create mode 100644 examples/opensees/respSpectra.py create mode 100644 examples/opensees/schematic.png create mode 100644 examples/opensees/short.tex create mode 100644 examples/opensees/velocity.input diff --git a/examples/opensees/DS_GenFunctionsV3.py b/examples/opensees/DS_GenFunctionsV3.py new file mode 100644 index 0000000..17dfe76 --- /dev/null +++ b/examples/opensees/DS_GenFunctionsV3.py @@ -0,0 +1,93 @@ +def DS_GetDir(cur_dir, t): + +# cur_dir = os.getcwd() + if "jupyter/MyData" in cur_dir: + cur_dir = cur_dir.split("MyData").pop() + input_dir = t.username + cur_dir + input_uri = "tapis://{}/{}".format(storage_id, input_dir) + input_uri = input_uri.replace(" ", "%20") + elif('jupyter/mydata' in cur_dir ): + cur_dir = cur_dir.split("myData").pop() + input_dir = t.username + cur_dir + input_uri = "tapis://{}/{}".format(storage_id, input_dir) + input_uri = input_uri.replace(" ", "%20") + + elif('jupyter/MyProjects' in cur_dir): + cur_dir = cur_dir.split("MyProjects/").pop() + PRJ = cur_dir.split("/")[0] + cur_dir = cur_dir.split(PRJ).pop() + import requests + + resp = requests.get( + f"https://designsafe-ci.org/api/projects/v2/{PRJ}", + headers={"x-tapis-token": t.access_token.access_token}, + ) + project_uuid = resp.json()["baseProject"]["uuid"] + input_dir = cur_dir + input_uri = "tapis://project-{}{}".format(project_uuid, cur_dir) + input_uri = input_uri.replace(" ", "%20") + elif "jupyter/projects" in cur_dir: + cur_dir = cur_dir.split("projects/").pop() + PRJ = cur_dir.split("/")[0] + cur_dir = cur_dir.split(PRJ).pop() + import requests + + resp = requests.get( + f"https://designsafe-ci.org/api/projects/v2/{PRJ}", + headers={"x-tapis-token": t.access_token.access_token}, + ) + project_uuid = resp.json()["baseProject"]["uuid"] + input_dir = cur_dir + input_uri = "tapis://project-{}{}".format(project_uuid, cur_dir) + input_uri = input_uri.replace(" ", "%20") + elif "jupyter/CommunityData" in cur_dir: + cur_dir = cur_dir.split("jupyter/CommunityData").pop() + input_dir = cur_dir + input_uri = "tapis://designsafe.storage.community/{}".format(input_dir) + input_uri = input_uri.replace(" ", "%20") + + return input_uri + +def DS_GetStatus(t, mjobUuid, tlapse = 15): + + import time + print(" Job launched. Status provided below") + print( + " Can also check in DesignSafe portal under - Workspace > Tools & Application > Job Status" + ) + + status = t.jobs.getJobStatus(jobUuid=mjobUuid).status + previous = "" + while True: + if status in ["FINISHED","FAILED","STOPPED"]: + break + status = t.jobs.getJobStatus(jobUuid=mjobUuid).status + if status == previous: + continue + else : + previous = status + print(f"\tStatus: {status}") + time.sleep(tlapse) + return status + +def DS_GetRuntime(t, mjobUuid): + + from datetime import datetime + + print("\nRuntime Summary") + print("---------------") + hist = t.jobs.getJobHistory(jobUuid=mjobUuid) + + time1 = datetime.strptime(hist[-1].created, "%Y-%m-%dT%H:%M:%S.%fZ") + time0 = datetime.strptime(hist[0].created, "%Y-%m-%dT%H:%M:%S.%fZ") + print("TOTAL time:", time1 - time0) + + for i in range(len(hist)): + if hist[i].eventDetail == 'RUNNING' : + time1 = datetime.strptime(hist[i+1].created, "%Y-%m-%dT%H:%M:%S.%fZ") + time0 = datetime.strptime(hist[i].created, "%Y-%m-%dT%H:%M:%S.%fZ") + print("RUNNING time:", time1 - time0) + if hist[i].eventDetail == 'QUEUED' : + time1 = datetime.strptime(hist[i+1].created, "%Y-%m-%dT%H:%M:%S.%fZ") + time0 = datetime.strptime(hist[i].created, "%Y-%m-%dT%H:%M:%S.%fZ") + print("QUEUED time:", time1 - time0) \ No newline at end of file diff --git a/examples/opensees/N10_T3.tcl b/examples/opensees/N10_T3.tcl new file mode 100644 index 0000000..5fd8a9b --- /dev/null +++ b/examples/opensees/N10_T3.tcl @@ -0,0 +1,515 @@ +# ######################################################### +# # +# Effective stress site response analysis for a layered # +# soil profile located on a level ground and underlain by # +# an elastic half-space. SSPquadUP elements are used. # +# The finite rigidity of the elastic half space is # +# considered through the use of a viscous damper at the # +# base. PM4Sand model (Bounlanger and Ziotopoulou 2017) # +# is used to simualte the liquefiable layer. # +# # +# Created by: Chris McGann # +# Pedro Arduino # +# Modified by: Long Chen # +# Andrew Makdisi # +# Alborz Ghofrani # +# --University of Washington-- # +# # +# ---> Basic units are kN and m (unless specified) # +# # +# ######################################################### + +wipe + +# --------------------------------------------------------- +# 1. DEFINE SOIL AND MESH GEOMETRY +# --------------------------------------------------------- + +# ---SOIL GEOMETRY +# number of soil layers +set numLayers 3 +# layer thicknessess (m) +set layerThick(3) 2.0 +set layerThick(2) 3.0 +set layerThick(1) 1.0 + +# depth of water table, dry +set waterTable 2.0 + +# ---MESH GEOMETRY +# number of elements in horizontal direction +set nElemX 1 +# number of nodes in horizontal direction +set nNodeX [expr $nElemX + 1] +# horizontal element size (m) +set sElemX 0.50 +# number of elements in vertical direction for each layer +set nElemY(3) 4 +set nElemY(2) 6 +set nElemY(1) 2 + +# define grade of slope (%) +set grade 0.0 +set g -9.81 + +set N160 10.0 +set Cd 46.0 +set Dr [expr {sqrt($N160 / $Cd)} ] +set Gs 2.67 +set emax 0.8 +set emin 0.5 +set void [expr $emax - $Dr * ($emax - $emin)] +set por [expr $void / (1 + $void)] +set rho_d [expr $Gs / (1 + $void)] +set rho_s [expr $rho_d *(1.0+$void/$Gs)] + +set K0 0.5 +set nu [expr $K0 / (1 + $K0)] + +# define properties of the underlying rock +set rockVS 182.0 +set rockDen $rho_s + +# ---GROUND MOTION PARAMETERS +# define velocity time history file +set velocityFile velocity.input +# time step in ground motion record +set motionDT 0.005 +# number of steps in ground motion record +set motionSteps 7998 + +# ---RAYLEIGH DAMPING PARAMETERS +set pi 3.141592654 +set fmin 5.01 +set Omegamin [expr $fmin * 2.0 * $pi] +set ximin 0.025 + +# factor to mass matrix +set a0 [expr $ximin * $Omegamin] +# factor to stiffness matrix +set a1 [expr $ximin / $Omegamin] + + +# calculate the thickness of soil profile +set soilThick 0.0 +for {set i 1} {$i <= $numLayers} {incr i} { + set soilThick [expr $soilThick + $layerThick($i)] +} + +# total number of elements in vertical direction +set nElemT 0 +set layerBound(1) 0 +for {set i 1} {$i <= $numLayers} {incr i} { + incr nElemT [expr $nElemY($i)*$nElemX] + set sElemY($i) [expr $layerThick($i)/$nElemY($i)] + set layerBound([expr $i+1]) [expr $layerBound($i) + $nElemY($i)] + puts "size: $sElemY($i)" + puts "layerBound $layerBound([expr $i+1])" +} +set layerBound(1) 1 + +# number of nodes in vertical direction in each layer +set nNodeT 0 +for {set k 1} {$k < $numLayers} {incr k 1} { + set nNodeL($k) [expr $nNodeX*$nElemY($k)] + puts "number of nodes in layer $k: $nNodeL($k)" + set nNodeT [expr $nNodeT + $nNodeL($k)] +} +set nNodeL($numLayers) [expr $nNodeX*($nElemY($numLayers) + 1)] +puts "number of nodes in layer $numLayers: $nNodeL($numLayers)" +set nNodeT [expr $nNodeT + $nNodeL($numLayers)] +puts "total number of nodes: $nNodeT" + +#----------------------------------------------------------------------------------------- +# 2. CREATE PORE PRESSURE NODES AND FIXITIES +#----------------------------------------------------------------------------------------- +model BasicBuilder -ndm 2 -ndf 3 + +set yCoord 0.0 +set count 0 +set gwt 1 +set waterHeight [expr $soilThick-$waterTable] +set nodesInfo [open nodesInfo.dat w] +# loop over layers +for {set k 1} {$k <= $numLayers} {incr k 1} { + # loop over nodes + for {set j 1} {$j <= $nNodeL($k)} {incr j $nNodeX} { + for {set i 1} {$i <= $nNodeX} {incr i} { + node [expr $j+$count+$i-1] [expr ($i-1)*$sElemX] $yCoord + puts $nodesInfo "[expr $j+$count+$i-1] [expr ($i-1)*$sElemX] $yCoord" + + # designate nodes above water table + if {$yCoord>=$waterHeight} { + set dryNode($gwt) [expr $j+$count+$i-1] + set gwt [expr $gwt+1] + } + } + + set yCoord [expr $yCoord + $sElemY($k)] + } + set count [expr $count + $nNodeL($k)] +} +close $nodesInfo + + +# define fixities for pore pressure nodes at base of soil column +for {set i 1} {$i <= $nNodeX} {incr i} { + fix $i 0 1 0 + # puts "fix $i 0 1 0" + if {$i > 1} { + equalDOF 1 $i 1 + # puts "equalDOF 1 $i 1" + } +} +puts "Finished creating all -ndf 3 boundary conditions..." + + +# define equal degrees of freedom for pore pressure nodes +for {set j [expr $nNodeX + 1]} {$j < $nNodeT} {incr j $nNodeX} { + for {set i $j} {$i < [expr $j + $nNodeX-1]} {incr i} { + equalDOF $j [expr $i+1] 1 2 + # puts "equalDOF $j [expr $i+1] 1 2" + } +} +puts "Finished creating equalDOF for pore pressure nodes..." + +# define pore pressure boundaries for nodes above water table +for {set i 1} {$i < $gwt} {incr i 1} { + fix $dryNode($i) 0 0 1 +} + +#----------------------------------------------------------------------------------------- +# 3. CREATE SOIL MATERIALS +#----------------------------------------------------------------------------------------- +set slope [expr atan($grade/100.0)] + +nDMaterial PM4Sand 3 $Dr 468.3 0.463 $rho_d 101.3 -1.00 $emax $emin 0.5 0.1 -1.0 -1.0 250.0 -1.00 33.0 $nu +puts "nDMaterial PM4Sand 2 $Dr 468.3 0.463 $rho_d 101.3 -1.00 $emax $emin 0.5 0.1 -1.0 -1.0 250.0 -1.00 33.0 $nu" +set thick(3) 1.0 +set xWgt(3) [expr $g*sin($slope)] +set yWgt(3) [expr $g*cos($slope)] +set uBulk(3) 2.2e6 +set hPerm(3) 1.0e-5 +set vPerm(3) 1.0e-5 +set eInit(3) $void + +nDMaterial PM4Sand 2 $Dr 584.1 0.450 $rho_s 101.3 -1.00 $emax $emin 0.5 0.1 -1.0 -1.0 250.0 -1.00 33.0 $nu +puts "nDMaterial PM4Sand 1 $Dr 584.1 0.450 $rho_s 101.3 -1.00 $emax $emin 0.5 0.1 -1.0 -1.0 250.0 -1.00 33.0 $nu" +set thick(2) 1.0 +set xWgt(2) [expr $g*sin($slope)] +set yWgt(2) [expr $g*cos($slope)] +set uBulk(2) 2.2e6 +set hPerm(2) 1.0e-5 +set vPerm(2) 1.0e-5 +set eInit(2) $void + +set E [expr 2 * $rho_s * $rockVS * $rockVS * (1 + 0.3)] +nDMaterial ElasticIsotropic 1 $E 0.3 $rho_s +set thick(1) 1.0 +set xWgt(1) [expr $g*sin($slope)] +set yWgt(1) [expr $g*cos($slope)] +set uBulk(1) 2.2e6 +set hPerm(1) 1.0e-9 +set vPerm(1) 1.0e-9 +set eInit(1) [expr 0.2 / (1 - 0.2)] + +puts "Finished creating all soil materials..." + +#----------------------------------------------------------------------------------------- +# 4. CREATE SOIL ELEMENTS +#----------------------------------------------------------------------------------------- +set elemInfo [open elementInfo.dat w] +set count 0 +for {set i 1} {$i <= $numLayers} {incr i 1} { + for {set j 1} {$j <= $nElemY($i)} {incr j 1} { + for {set k 1} {$k <= $nElemX} {incr k} { + set nI [expr ($nNodeX)*($j+$count-1) + $k] + set nJ [expr $nI + 1] + set nK [expr $nI + $nNodeX + 1] + set nL [expr $nI + $nNodeX] + + # permeabilities are initially set at 1.0 m/s for gravity analysis, values are updated post-gravity + element SSPquadUP [expr ($nElemX)*($j+$count-1) + $k] $nI $nJ $nK $nL $i $thick($i) $uBulk($i) 1.0 1.0 1.0 $eInit($i) 1.0e-8 $xWgt($i) $yWgt($i) + puts $elemInfo "[expr ($nElemX)*($j+$count-1) + $k] $nI $nJ $nK $nL $i" + } + } + set count [expr $count + $nElemY($i)] +} +close $elemInfo +puts "Finished creating all soil elements..." + +#----------------------------------------------------------------------------------------- +# 6. LYSMER DASHPOT +#----------------------------------------------------------------------------------------- +model BasicBuilder -ndm 2 -ndf 2 + +# define dashpot nodes +set dashF [expr $nNodeT+1] +set dashS [expr $nNodeT+2] + +node $dashF 0.0 0.0 +node $dashS 0.0 0.0 + +# define fixities for dashpot nodes +fix $dashF 1 1 +fix $dashS 0 1 + +# define equal DOF for dashpot and base soil node +equalDOF 1 $dashS 1 +puts "Finished creating dashpot nodes and boundary conditions..." + +# define dashpot material +set colArea [expr $sElemX*$thick(1)] +set dashpotCoeff [expr $rockVS*$rockDen] +uniaxialMaterial Viscous [expr $numLayers+1] [expr $dashpotCoeff*$colArea] 1 + +# define dashpot element +element zeroLength [expr $nElemT+1] $dashF $dashS -mat [expr $numLayers+1] -dir 1 +puts "Finished creating dashpot material and element..." + +#----------------------------------------------------------------------------------------- +# 7. CREATE GRAVITY RECORDERS +#----------------------------------------------------------------------------------------- + +# create list for pore pressure nodes +set nodeList3 {} +set channel [open "nodesInfo.dat" r] +set count 0; +foreach line [split [read -nonewline $channel] \n] { +set count [expr $count+1]; +set lineData($count) $line +set nodeNumber [lindex $lineData($count) 0] +lappend nodeList3 $nodeNumber +} +close $channel + +# record nodal displacment, acceleration, and porepressure +eval "recorder Node -file Gdisplacement.out -time -node $nodeList3 -dof 1 2 disp" +eval "recorder Node -file Gacceleration.out -time -node $nodeList3 -dof 1 2 accel" +eval "recorder Node -file GporePressure.out -time -node $nodeList3 -dof 3 vel" +# record elemental stress and strain (files are names to reflect GiD gp numbering) +recorder Element -file Gstress.out -time -eleRange 1 $nElemT stress 3 +recorder Element -file Gstrain.out -time -eleRange 1 $nElemT strain +puts "Finished creating gravity recorders..." + + +# ----------------------------------------------------------------------------------------- +# 8. CREATE GID FLAVIA.MSH FILE FOR POSTPROCESSING +# ----------------------------------------------------------------------------------------- + +set meshFile [open freeFieldEffective.flavia.msh w] +puts $meshFile "MESH ffBrick dimension 2 ElemType Quadrilateral Nnode 4" +puts $meshFile "Coordinates" +puts $meshFile "#node_number coord_x coord_y" +set yCoord 0.0 +set count 0 +# loop over layers +for {set k 1} {$k <= $numLayers} {incr k 1} { + # loop over nodes + for {set j 1} {$j <= $nNodeL($k)} {incr j $nNodeX} { + for {set i 1} {$i <= $nNodeX} {incr i} { + puts $meshFile "[expr $j+$count+$i-1] [expr ($i-1)*$sElemX] $yCoord" + } + set yCoord [expr $yCoord + $sElemY($k)] + } + set count [expr $count + $nNodeL($k)] +} +puts $meshFile "end coordinates" +puts $meshFile "Elements" +puts $meshFile "# element node1 node2 node3 node4" +set count 0 +# loop over layers +for {set i 1} {$i <= $numLayers} {incr i 1} { + for {set j 1} {$j <= $nElemY($i)} {incr j 1} { + for {set k 1} {$k <= $nElemX} {incr k} { + set nI [expr ($nNodeX)*($j+$count-1) + $k] + set nJ [expr $nI + 1] + set nK [expr $nI + $nNodeX + 1] + set nL [expr $nI + $nNodeX] + puts $meshFile "[expr $j+$count] $nI $nJ $nK $nL" + } + } + set count [expr $count + $nElemY($i)] +} +puts $meshFile "end elements" +close $meshFile + + +#----------------------------------------------------------------------------------------- +# 9. DEFINE ANALYSIS PARAMETERS +#----------------------------------------------------------------------------------------- + +#---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION + +# duration of ground motion (s) +set duration [expr $motionDT*$motionSteps] + +set AnalysisdT 0.001 +set nSteps [expr int(floor($duration/$AnalysisdT)+1)] +set dT $AnalysisdT + +puts "number of steps in analysis: $nSteps" +puts "analysis time step: $dT" + +#---ANALYSIS PARAMETERS +# Newmark parameters + +set gamma [expr 5.0/6.0] +set beta [expr 4.0/9.0] + +#----------------------------------------------------------------------------------------- +# 10. GRAVITY ANALYSIS +#----------------------------------------------------------------------------------------- + +# update materials to ensure elastic behavior +for {set i 1} {$i <= $numLayers} {incr i} { + updateMaterialStage -material $i -stage 0 +} + +# fix bottom nodes for additional stability +# fix 1 1 1 0 + +constraints Transformation +test NormDispIncr 1e-4 35 1 +algorithm Newton +numberer RCM +#system SparseGeneral +system ProfileSPD +integrator Newmark $gamma $beta +analysis Transient + +set startT [clock seconds] +analyze 10 1.0 +puts "Finished with elastic gravity analysis..." + +# update materials to consider plastic behavior +for {set i 1} {$i <= $numLayers} {incr i} { + updateMaterialStage -material $i -stage 1 +} +for {set i 1} {$i <= $numLayers} {incr i} { + setParameter -value 0 -eleRange $layerBound($i) $layerBound([expr $i+1]) FirstCall $i + setParameter -value 0.3 -eleRange $layerBound($i) $layerBound([expr $i+1]) poissonRatio $i +} +# plastic gravity loading +analyze 10 1.0 + +# remove extra bottom fixity for dynamic analysis +# remove sp 1 1 +puts "Finished with plastic gravity analysis..." + +#----------------------------------------------------------------------------------------- +# 11. UPDATE ELEMENT PERMEABILITY VALUES FOR POST-GRAVITY ANALYSIS +#----------------------------------------------------------------------------------------- + +# update permeability parameters for each element +for {set i 1} {$i <= $numLayers} {incr i} { + setParameter -value [expr $hPerm($i) / 9.81] -eleRange $layerBound($i) $layerBound([expr $i+1]) hPerm + setParameter -value [expr $vPerm($i) / 9.81] -eleRange $layerBound($i) $layerBound([expr $i+1]) vPerm +} +puts "Finished updating permeabilities for dynamic analysis..." + +#----------------------------------------------------------------------------------------- +# 12. CREATE POST-GRAVITY RECORDERS +#----------------------------------------------------------------------------------------- + +# reset time and analysis +setTime 0.0 +wipeAnalysis +remove recorders + +# recorder time step +set recDT 0.01 + +# record nodal displacment, acceleration, and porepressure +eval "recorder Node -file displacement.out -time -dT $recDT -node $nodeList3 -dof 1 2 disp" +eval "recorder Node -file acceleration.out -time -dT $recDT -node $nodeList3 -dof 1 2 accel" +eval "recorder Node -file porePressure.out -time -dT $recDT -node $nodeList3 -dof 3 vel" +# record elemental stress and strain (files are names to reflect GiD gp numbering) +recorder Element -file stress.out -time -dT $recDT -eleRange 1 $nElemT stress 3 +recorder Element -file strain.out -time -dT $recDT -eleRange 1 $nElemT strain +puts "Finished creating all recorders..." + +#----------------------------------------------------------------------------------------- +# 13. DYNAMIC ANALYSIS +#----------------------------------------------------------------------------------------- + +model BasicBuilder -ndm 2 -ndf 3 + +# define constant scaling factor for applied velocity +set cFactor [expr $colArea*$dashpotCoeff] + +# timeseries object for force history +set mSeries "Path -dt $motionDT -filePath $velocityFile -factor $cFactor" + +# loading object +pattern Plain 10 $mSeries { + load 1 1.0 0.0 0.0 +} +puts "Dynamic loading created..." + +set gamma 0.5 +set beta 0.25 + + +constraints Transformation +test NormDispIncr 1.0e-4 35 1 +algorithm Newton +numberer RCM +#system SparseGeneral +system ProfileSPD +integrator Newmark $gamma $beta +rayleigh $a0 $a1 0.0 0.0 +analysis Transient + +# Analyze and use substepping if needed +set remStep $nSteps +set success 0 + +proc subStepAnalyze {dT subStep} { + if {$subStep > 10} { + return -10 + } + for {set i 1} {$i < 3} {incr i} { + puts "Try dT = $dT" + set success [analyze 1 $dT] + if {$success != 0} { + set success [subStepAnalyze [expr $dT/2.0] [expr $subStep+1]] + if {$success == -10} { + puts "Did not converge." + return $success + } + } else { + if {$i==1} { + puts "Substep $subStep : Left side converged with dT = $dT" + } else { + puts "Substep $subStep : Right side converged with dT = $dT" + } + } + } + return $success +} + +puts "Start analysis" +set startT [clock seconds] + +while {$success != -10} { + set subStep 0 + set success [analyze $remStep $dT] + if {$success == 0} { + puts "Analysis Finished" + break + } else { + set curTime [getTime] + puts "Analysis failed at $curTime . Try substepping." + set success [subStepAnalyze [expr $dT/2.0] [incr subStep]] + set curStep [expr int($curTime/$dT + 1)] + set remStep [expr int($nSteps-$curStep)] + puts "Current step: $curStep , Remaining steps: $remStep" + } +} +set endT [clock seconds] +puts "loading analysis execution time: [expr $endT-$startT] seconds." + +puts "Finished with dynamic analysis..." + +wipe diff --git a/examples/opensees/ShortReport.tex b/examples/opensees/ShortReport.tex new file mode 100644 index 0000000..29fe5d3 --- /dev/null +++ b/examples/opensees/ShortReport.tex @@ -0,0 +1,98 @@ +%######################################################## +% +% Postprocessing LaTeX script +% +% Copyright: UW Computational Mechanics Group +% Pedro Arduino +% +% Participants: Alborz Ghofrani +% Long Chen +% +%------------------------------------------------------- + +\documentclass[11pt,fleqn]{article} + +\usepackage[T1]{fontenc} +\usepackage[ansinew]{inputenc} +\usepackage{fullpage, url} +\usepackage{amsmath, amsfonts} + +\usepackage[normalem]{ulem} +\usepackage{longtable} +\usepackage{xcolor} +\usepackage{graphicx} +\newcommand\suppress[1]{} +\newcommand\deleted[1]{\xout{#1}} +\newcommand\revised[1]{\uline{#1}} +\newlength\wvtextpercent +\setlength\wvtextpercent{0.009\textwidth} + +% \newif\ifpdf +% \ifx\pdfoutput\undefined +% \pdffalse +% \else +% \pdfoutput=1 +% \pdftrue +% \fi + +%\ifpdf +% \usepackage[pdftex]{xcolor} +% \usepackage[pdftex]{graphicx} +% \pdfinfo { +% /Title (Materials Modeling) +% /Subject (Transition from 1D to 3D) +% /Author (Peter Mackenzie) +% /Keywords (CEE503) +% } +%\else +% \usepackage{xcolor} +% \usepackage{graphicx} +%\fi + +\setlength{\parindent}{0pt} +\setlength{\parskip}{1ex plus 0.5ex minus 0.2ex} + +\input{short} +\input{macros} + + +\begin{document} + +%\tableofcontents +%\newpage + \begin{center} + + \textbf{{DesignSafe Example\hfill{}1D Site Response Examples\hfill{}May 2017}} + + \end{center} + + +\section{Soil profile} + +Free field response of single soil profile subjected to earthquake excitation + +\begin{figure}[h!] +\centering +\includegraphics[scale=0.7]{schematic} +\caption{1D Model} +\end{figure} + +\newpage +\section{Results} + +% Plot Results +\begin{figure}[h!] +\centering +\includegraphics[width=5in]{surfaceAccel} +\caption{Surface acceleration time history} +\end{figure} + +\begin{figure}[h!] +\centering +\includegraphics[width=5in]{logSpectra} +\caption{Surface response spectra} +\end{figure} + + + +\end{document} diff --git a/examples/opensees/interactiveplot.py b/examples/opensees/interactiveplot.py new file mode 100644 index 0000000..8167090 --- /dev/null +++ b/examples/opensees/interactiveplot.py @@ -0,0 +1,73 @@ +""" +Create interactive plot for disp and pwp +@author: Long Chen +""" + +import matplotlib.pyplot as plt +from matplotlib import gridspec +from ipywidgets import interactive +import ipywidgets as widgets +import numpy as np + +def pwpplot(timeStep): + Step = int(timeStep / 0.01)-1 + plt.subplot(211) + plt.plot(time, uu) + plt.plot(time[Step],uu[Step],'ro') + plt.ylabel('pwp(kPa)') + plt.grid() + plt.subplot(212) + plt.plot(time,acc_input) + plt.plot(time[Step],acc_input[Step],'ro') + plt.xlabel('time(s)') + plt.ylabel('acceleration(g)') + plt.grid() + +def dispplot(timeStep): + Step = int(timeStep / 0.01)-1 + plt.figure(figsize=(7, 8)) + ax0 = plt.subplot(gs[0]) + ax0.plot(maxdisp[0, ::2], nodes[::2, 2], 'b--') + ax0.plot(mindisp[0, ::2], nodes[::2, 2], 'b--') + ax0.plot(disp[Step, ::4], nodes[::2, 2]) + plt.xlabel('displacement(m)') + plt.ylabel('Elevation(m)') + plt.grid() + ax1 = plt.subplot(gs[1]) + ax1.plot(time,acc_input) + ax1.plot(time[Step],acc_input[Step],'ro') + plt.xlabel('time(s)') + plt.ylabel('acceleration(g)') + plt.grid() + +def createpwpplot(): + global time, acc_input, uu + pwp = np.loadtxt('porePressure.out') + time = pwp[:,0] + pwp = np.delete(pwp, 0, 1) + uexcess = pwp - pwp[0, :] + uu = uexcess[0:len(time), 12] + acc = np.loadtxt('acceleration.out') + acc_input = acc[:, 1] + + return interactive(pwpplot,timeStep = widgets.FloatSlider(min = 0.01, max = time[-1], step = 0.01)) + + +def createDispplot(): + global maxdisp, mindisp, nodes, disp, gs + nodes = np.loadtxt('nodesInfo.dat') + disp = np.loadtxt('displacement.out') + disp = np.delete(disp, 0, 1) + disp = (disp.transpose() - disp[:,0]).transpose() + ndof = 2 + nnodes = nodes.shape[0] + maxdisp = np.amax(disp, axis=0) + mindisp = np.amin(disp, axis=0) + maxdisp = maxdisp.reshape(ndof, nnodes, order="F") + mindisp = mindisp.reshape(ndof, nnodes, order="F") + gs = gridspec.GridSpec(2, 1, height_ratios=[6, 1]) + + return interactive(dispplot,timeStep = widgets.FloatSlider(min = 0.01, max = time[-1], step = 0.01), continuous_update=False) + +if __name__ == "__main__": + createpwpplot() diff --git a/examples/opensees/opensees-dapi-v3.ipynb b/examples/opensees/opensees-dapi-v3.ipynb new file mode 100644 index 0000000..8baa4c7 --- /dev/null +++ b/examples/opensees/opensees-dapi-v3.ipynb @@ -0,0 +1,577 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Free Field Analysis Example\n", + "Pedro Arduino - UW Computational Geomechanics Group" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example shows how to run OpenSees in DesignSafe from a jupyter notebook and how to postprocess the output results using python scripts, generate a LaTex report, and create interactive plots. \n", + "\n", + "A simple 1D free field analysis of a liquefiable soil layer is analyzed using OpenSees. An schematic of the soil profile in shown in the Figure. The soil profile consists of a 1 m dry crust, 3 m liquefiable layer, and 1 m of elastic base. The ground water table is at 2 m. An earthquake excitation is applied at the bottom of the soil column. A compliant rock is considered in the analysis. \n", + "\n", + "The results are presented in terms of:\n", + "\n", + "a) Time history of acceleration at the surface and corresponding response spectra.\n", + "\n", + "b) Profiles of maximum displacement, peak horizontal acceleration (PHA), maximum shear strain, and stress ratio\n", + "\n", + "c) Stress strain plots for a point near the center of the liquefiable zone, and\n", + "\n", + "d) Evolution of pore water pressure for a point near the center of the liquefiable zone. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Setup agave and start OpenSees job" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Setup job description" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "!pip install --user --upgrade setuptools --quiet\n", + "!pip install --user --only-binary=:all: atomicwrites==1.4.0 --quiet\n", + "!pip install --user \"jsonschema<4.18.0\" --quiet\n", + "!pip install git+https://github.com/DesignSafe-CI/dapi.git@tapisv3 --user --quiet" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import dapi\n", + "import uuid\n", + "from datetime import date\n", + "import json" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Authenticate\n", + "t = dapi.auth.init()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Define inputs\n", + "cur_dir = os.getcwd()\n", + "input_uri = dapi.jobs.get_ds_path_uri(t, cur_dir)\n", + "input_filename = \"N10_T3.tcl\"" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "---Job Info---\n", + "\n", + "{\n", + " \"name\": \"opensees-express_20241001_174715\",\n", + " \"appId\": \"opensees-express\",\n", + " \"appVersion\": \"3.7.0\",\n", + " \"execSystemId\": \"wma-exec-01\",\n", + " \"maxMinutes\": 30,\n", + " \"archiveOnAppError\": true,\n", + " \"fileInputs\": [\n", + " {\n", + " \"name\": \"Input Directory\",\n", + " \"sourceUrl\": \"tapis://designsafe.storage.default/kks32/freeFieldEffectiveJupyter\"\n", + " }\n", + " ],\n", + " \"execSystemLogicalQueue\": \"development\",\n", + " \"nodeCount\": 1,\n", + " \"coresPerNode\": 1,\n", + " \"parameterSet\": {\n", + " \"envVariables\": [\n", + " {\n", + " \"key\": \"tclScript\",\n", + " \"value\": \"N10_T3.tcl\"\n", + " }\n", + " ]\n", + " }\n", + "}\n" + ] + } + ], + "source": [ + "job_info = dapi.jobs.generate_job_info(t, 'opensees-express', input_uri, input_filename)\n", + "job_info['maxMinutes'] = 30\n", + "job_info['execSystemLogicalQueue'] = 'development'\n", + "print(\"\\n---Job Info---\\n\\n\" + json.dumps(job_info, indent=2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Run job" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "# Submit job\n", + "job = t.jobs.submitJob(**job_info)\n", + "jobUuid=job.uuid" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Waiting for job to start: 3it [00:45, 15.07s/it, Status: RUNNING] \n", + "Monitoring job: 0%| | 0/120 [00:00 Fix\n", + "\n", + "```json\n", + "job_description[\"archiveSystemId\"] = \"designsafe.storage.default\"\n", + "job_description[\"archiveSystemDir\"] = (\n", + " f\"{t.username}/tapis-jobs-archive/${{JobCreateDate}}/${{JobName}}-${{JobUUID}}\"\n", + ")\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "# %cd ..\n", + "cur_dir_name = cur_dir.split('/').pop() \n", + "os.chdir(jobinfo.archiveSystemDir.replace(user,'/home/jupyter/MyData'))\n", + "if not os.path.exists(cur_dir_name):\n", + " os.makedirs(cur_dir_name)\n", + "os.chdir(cur_dir_name) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import python libraries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot acceleration time history" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot acceleration time hisotory and response spectra on log-linear scale" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from plotAcc import plot_acc\n", + "plot_acc()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot profiles" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot profiles of max displacement, PGA, max shear strain, stress ratio and plot stress strain near the center of liquefiable layer " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "from plotProfile import plot_profile\n", + "plot_profile()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot excess pore water pressure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from plotPorepressure import plot_porepressure\n", + "plot_porepressure()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Generate LaTeX Report " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#os.system('/usr/bin/pdflatex -interaction nonstopmode ShortReport.tex')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Before we start let us install a python package for plotting\n", + "try:\n", + " import rst2pdf\n", + "\n", + "except:\n", + " import pip\n", + " pip.main(['install', '--user', 'rst2pdf'])\n", + " print(\"********* Please restart the session ***********\")\n", + " \n", + "import rst2pdf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "!{sys.executable} -m pip install rst2pdf -q" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 2024 - JupyterHub\n", + "os.system('/home/jupyter/.local/bin/rst2pdf ShortReport.rst ShortReport.pdf')\n", + "# 2022 - JupyterHub\n", + "# os.system('/opt/conda/bin/rst2pdf ShortReport.rst ShortReport.pdf')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class PDF(object):\n", + " def __init__(self, pdf, size=(200,200)):\n", + " self.pdf = pdf\n", + " self.size = size\n", + "\n", + " def _repr_html_(self):\n", + " return ''.format(self.pdf, self.size)\n", + "\n", + " def _repr_latex_(self):\n", + " return r'\\includegraphics[width=1.0\\textwidth]{{{0}}}'.format(self.pdf)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# pdf_fn = jobinfo.archiveSystemDir.replace(user, '/user/' + user + '/files/MyData')\n", + "pdf_fn = jobinfo.archiveSystemDir.replace('/'+user, '../../../MyData')\n", + "\n", + "pdf_fn += '/'\n", + "pdf_fn += cur_dir.split('/')[-1]\n", + "pdf_fn += '/ShortReport.pdf'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "PDF(pdf_fn , (750,600))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Create Interactive Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pore water pressure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from interactiveplot import createpwpplot, createDispplot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "createpwpplot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Displacement" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "createDispplot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "IMAGE_NAME": "taccsciapps/ds-nb-img:base-0.1.9", + "UUID": "73e0880d-9b87-11ec-9c1c-13579dd95994", + "celltoolbar": "Raw Cell Format", + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/opensees/plotAcc.py b/examples/opensees/plotAcc.py new file mode 100644 index 0000000..d812f3c --- /dev/null +++ b/examples/opensees/plotAcc.py @@ -0,0 +1,64 @@ +######################################################### +# +# Postprocessing python script +# +# Copyright: UW Computational Mechanics Group +# Pedro Arduino +# +# Participants: Alborz Ghofrani +# Long Chen +# +#------------------------------------------------------- + +import numpy as np +import matplotlib.pyplot as plt + +from respSpectra import resp_spectra + + +def plot_acc(ndof=2): + """ + Plot acceleration time history and response spectra + """ + acc = np.loadtxt('acceleration.out') + time = acc[:, 0] + acc = np.delete(acc, 0, 1) + # Bandaid to remove last 2 nodes associated with dashpot (not for all Openees) + #acc = acc[:,0:-4] + + [nstep, temp] = acc.shape + nnode = int(temp/ndof) + a = acc.reshape(nstep, ndof, nnode, order="F") / 9.81 + + # plot acceleration time history + plt.figure() + plt.plot(time, a[:, 0, nnode-1]) + plt.grid() + plt.xlabel('time (sec)') + plt.ylabel('acceleration (g)') + plt.savefig('surfaceAccel.eps') + plt.savefig('surfaceAccel.png') + + # build response spectra + [p, umax, vmax, amax] = resp_spectra(a[:, 0, nnode-1], time[-1], nstep) + + # response spectra on log-linear plot + plt.figure() + plt.subplot(3, 1, 1) + plt.semilogx(p, amax) + plt.grid() + plt.ylabel('$S_a (g)$') + plt.subplot(3, 1, 2) + plt.semilogx(p, vmax) + plt.grid() + plt.ylabel('$S_v (m/s)$') + plt.subplot(3, 1, 3) + plt.semilogx(p, umax) + plt.grid() + plt.ylabel('$S_d (m)$') + plt.xlabel('$Period (s)$') + plt.savefig('logSpectra.eps') + plt.savefig('logSpectra.png') + +if __name__ == "__main__": + plot_acc() diff --git a/examples/opensees/plotPorepressure.py b/examples/opensees/plotPorepressure.py new file mode 100644 index 0000000..c34c7eb --- /dev/null +++ b/examples/opensees/plotPorepressure.py @@ -0,0 +1,35 @@ +######################################################### +# +# Postprocessing python script +# +# Copyright: UW Computational Mechanics Group +# Pedro Arduino +# +# Participants: Alborz Ghofrani +# Long Chen +# +#------------------------------------------------------- + +import numpy as np +import matplotlib.pyplot as plt + + +def plot_porepressure(): + """ + Plot pore water pressure + """ + porepressure = np.loadtxt('porePressure.out') + time = porepressure[:, 0] + porepressure = np.delete(porepressure, 0, 1) + uexcess = porepressure - porepressure[0, :] + + plt.figure() + plt.plot(time, uexcess[:, 12]) + plt.xlabel('Time(s)') + plt.ylabel('u_excess(kPa)') + plt.grid() + plt.savefig('porePressure.eps') + plt.savefig('porePressure.png') + +if __name__ == "__main__": + plot_porepressure() diff --git a/examples/opensees/plotProfile.py b/examples/opensees/plotProfile.py new file mode 100644 index 0000000..2ab9798 --- /dev/null +++ b/examples/opensees/plotProfile.py @@ -0,0 +1,87 @@ +######################################################### +# +# Postprocessing python script +# +# Copyright: UW Computational Mechanics Group +# Pedro Arduino +# +# Participants: Alborz Ghofrani +# Long Chen +# +#------------------------------------------------------- + +from plotStressStrain import plot_stress_strain + +import numpy as np +import matplotlib.pyplot as plt + + +def plot_profile(ndof=2, nstraincomp=3, nstresscomp=3): + """ + Plot maximum displacement, PGA and maximum shear strain and maximum cyclic stress ratio + """ + nodes = np.loadtxt('nodesInfo.dat') + disp = np.loadtxt('displacement.out') + acc = np.loadtxt('acceleration.out') + strain = np.loadtxt('strain.out') + stress = np.loadtxt('stress.out') + + time = acc[:,0] + disp = np.delete(disp, 0, 1) + acc = np.delete(acc, 0, 1) + strain = np.delete(strain, 0, 1) + stress = np.delete(stress, 0, 1) + + # Bandaid to remove last 2 nodes associated with dashpot (not correct for all Openees) + #disp = disp[:,0:-4] + #acc = acc[:,0:-4] + + # subtact base displacement + disp = (disp.transpose() - disp[:,0]).transpose() + maxdisp = np.amax(np.abs(disp), axis=0) + pga = np.amax(np.abs(acc), axis=0) + maxstrain = np.amax(np.abs(strain), axis=0) + maxstressratio = np.amax(np.abs(stress[:, 2::nstresscomp]), axis=0) + maxstressratio = maxstressratio / np.abs(stress[0, 1::nstresscomp]) + + [nstep, temp] = strain.shape + nelem = int(temp / nstraincomp) + nnodes = nodes.shape[0] + + stress = stress.reshape(nstep, nstresscomp, nelem, order="F") + strain = strain.reshape(nstep, nstraincomp, nelem, order="F") + maxdisp = maxdisp.reshape(ndof, nnodes, order="F") + pga = pga.reshape(ndof, nnodes, order="F") + maxstrain = maxstrain.reshape(nstraincomp, nelem, order="F") + + plt.figure(figsize=(12, 6)) + plt.subplot(1, 4, 1) + plt.plot(maxdisp[0, ::2], nodes[::2, 2]) + plt.xticks(np.arange(0.0, max(maxdisp[0, ::2]), 0.2)) + plt.grid() + plt.xlabel('Maximum Displacement(m)') + plt.ylabel('Elevation(m)') + + plt.subplot(1, 4, 2) + plt.plot(pga[0, ::2] / 9.81, nodes[::2, 2]) + plt.xticks(np.arange(0.0, max(pga[0, ::2]) / 9.81, 0.2)) + plt.grid() + plt.xlabel('PHA(g)') + + plt.subplot(1, 4, 3) + plt.plot(maxstrain[2, :]*100.0, nodes[:-2:2, 2] + np.diff(nodes[::2, 2])) + plt.grid() + plt.xlabel('Maximum Shear Strain(%)') + + plt.subplot(1, 4, 4) + plt.plot(maxstressratio, nodes[:-2:2, 2] + np.diff(nodes[::2, 2])) + plt.xticks(np.arange(0.0, max(maxstressratio), 0.2)) + plt.grid() + plt.xlabel('$(\\tau / \sigma_{v0})_{max} $') + plt.savefig('profilePlot.eps') + plt.savefig('profilePlot.png') + + plot_stress_strain(stress[:, 2, 6], strain[:, 2, 6]) + +if __name__ == "__main__": + plot_profile() diff --git a/examples/opensees/plotStressStrain.py b/examples/opensees/plotStressStrain.py new file mode 100644 index 0000000..4b410d5 --- /dev/null +++ b/examples/opensees/plotStressStrain.py @@ -0,0 +1,26 @@ +######################################################### +# +# Postprocessing python script +# +# Copyright: UW Computational Mechanics Group +# Pedro Arduino +# +# Participants: Alborz Ghofrani +# Long Chen +# +#------------------------------------------------------- + +import matplotlib.pyplot as plt + + +def plot_stress_strain(stress,strain): + """ + Plot stress strain curve + """ + plt.figure() + plt.plot(strain,stress) + plt.xlabel('strain(%)') + plt.ylabel('stress(kPa)') + plt.grid() + plt.savefig('stressstrain.eps') + plt.savefig('stressstrain.png') diff --git a/examples/opensees/respSpectra.py b/examples/opensees/respSpectra.py new file mode 100644 index 0000000..991be80 --- /dev/null +++ b/examples/opensees/respSpectra.py @@ -0,0 +1,70 @@ +######################################################### +# +# Postprocessing python script +# +# Copyright: UW Computational Mechanics Group +# Pedro Arduino +# +# Participants: Alborz Ghofrani +# Long Chen +# +#------------------------------------------------------- + +import numpy as np + + +def resp_spectra(a, time, nstep): + """ + This function builds response spectra from acceleration time history, + a should be a numpy array,T and nStep should be integers. + """ + + # add initial zero value to acceleration and change units + a = np.insert(a, 0, 0) + # number of periods at which spectral values are to be computed + nperiod = 100 + # define range of considered periods by power of 10 + minpower = -3.0 + maxpower = 1.0 + # create vector of considered periods + p = np.logspace(minpower, maxpower, nperiod) + # incremental circular frequency + dw = 2.0 * np.pi / time + # vector of circular freq + w = np.arange(0, (nstep+1)*dw, dw) + # fast fourier Horm of acceleration + afft = np.fft.fft(a) + # arbitrary stiffness value + k = 1000.0 + # damping ratio + damp = 0.05 + umax = np.zeros(nperiod) + vmax = np.zeros(nperiod) + amax = np.zeros(nperiod) + # loop to compute spectral values at each period + for j in range(0, nperiod): + # compute mass and dashpot coeff to produce desired periods + m = ((p[j]/(2*np.pi))**2)*k + c = 2*damp*(k*m)**0.5 + h = np.zeros(nstep+2, dtype=complex) + # compute transfer function + for l in range(0, int(nstep/2+1)): + h[l] = 1./(-m*w[l]*w[l] + 1j*c*w[l] + k) + # mirror image of Her function + h[nstep+1-l] = np.conj(h[l]) + + # compute displacement in frequency domain using Her function + qfft = -m*afft + u = np.zeros(nstep+1, dtype=complex) + for l in range(0, nstep+1): + u[l] = h[l]*qfft[l] + + # compute displacement in time domain (ignore imaginary part) + utime = np.real(np.fft.ifft(u)) + + # spectral displacement, velocity, and acceleration + umax[j] = np.max(np.abs(utime)) + vmax[j] = (2*np.pi/p[j])*umax[j] + amax[j] = (2*np.pi/p[j])*vmax[j] + + return p, umax, vmax, amax diff --git a/examples/opensees/schematic.png b/examples/opensees/schematic.png new file mode 100644 index 0000000000000000000000000000000000000000..f8fd1353826cb11ce14e1a052beb9f4276b2b8dc GIT binary patch literal 7794 zcmd5>byQSupQi*wLO~j&K|#8^OF~3uV3c+khVDUX0Fe#>0R^N(fuROz>F%7NLAr(# zkdWB>vB?ANM@>dG5K-{odz~@24US^fk#y7)Wq%aLBZ^)C{rX8V(L_ zA~7EJ&R_Zxi5+k~3^i47%7&S?u@eGYWj$paoXQXPuB`50=XYJTKpr?aCNv-OJ{AQ$~Z1iGEXX=stfW&+EL$z#kRGT%e$H z%)z1uY{urGeZL*~WBedhw#09w-lFlW{=jFX-dKu;mX=I~Gd|~*3p`Qvgq?4g&68o* zU4cAYHK3T@-7g*vqJ3&!rWkwSB}WgvG+JmTD*lZqM#w-Z!^!tx>(?Jo*@uuTTgS)8 zqN1WBm_Pt+DGQmP{O|LRXSNUs1UzjkD z{)hx@&~6BasLzRi!lkXW4WN^C9o^oxYl?q5jsH83R*P4=3K5|Zp&g;?(8&iRNP&DsWIN9}Q4y!gjHFjx|tOb?46l1XdZ`;2SMDV_36|dX-QEM=K8lc@g7r0 zec~N6CH#&*9P{*RM@*ZNaVOy}YjszKB_oi!-xyn?m{*5ZAUIT)7+yirk4Fc&2kdoi zmAl%VRIhA<51sM^f3{~opp3%=23 ztPHb|-JGZzsQ2GFbuKT-*A?Wov{!9>AOoDDqJ5q(NKi=J-wKX&>)~ zxqC+r2wi<@NBQS}M<05()F6}t8P3tiO!GDSpbh6tIcD6UPqef-MZ=s$wv(G z4Gy;H(hY*e8qnY;AEjriuHt5)*sXp!WU-8SdJElvl7IpfwZX?Othy<{tuCSNFnvE$ zucfl5;{$AnjJG4=YU^}3g_#1}i*0CtavwI}V5oM2a(=$o$funyEDJ>?I?oEJFn9FhnvV18E$9RS0yeZ2kEg#GSbu?Wrx)XgQ)63MqRl*ibZx8=i&T<$R~B>l33OgY~Yt2mqc7WVsyI44T$BhUG`|N`&K3rQs6XY z54i6w8w1Q>`{CXiFbtJ}0lmAPrH@4$=Mne~2|x81WmSlZLwk$YUoxT;crv3EL$%1_ zy1LZm)5XxxFIGvlrlTCJ5+-K_%L0;e99l4k1%cJ;1o$zHV7Y@x!`#@d7i!-vwn=eC z)e59Hm0FO~d2w;k*Vm`SdYIh7bIoF|cs+dY&Ye5hRsF9o%M%_>?@dllc6N4#P*cMd z8HT{F^`-}86tpxn8|&+(NB+IV2vx9n_$$)%`#02bE&;uWlvU8RBp;R z<2S`(XTMPjvXJPrq@v<&&z;!+C0>>!f)_nY;*ZYx-#9zRAq~-DID0r( zem7+^!L3)%c?!Q478ddgQ7PQERTmvx*&h}7K`KQ;_-ynjM#&SXN#NJ}MJW^_!Y#UL z)$LO~Nb{=`qwaa#6B2LMByo8urns0GQvdxWiGzNLVe^5$;o)I2{Ak+3D!~!()%LJ7 zqHvUPsmK5=mObbplOaTTrFw*s6to$Zt#WsF&$cr+e*`hOslFY${M{r9qDf6nO-jc8 z23HB?(L_DHS4-1jDFtcgE3HWP!K>qf(0=Cb%^K+nyiP!CJ&HTgsX z*U)I`xUQCZ+ne(dTen|ZUe(Qt(I&h7uZ;M9Qt{z8t!G4W>=_$Hp1NP#4fe>;d)*D( zZW8fxXg)bsKT=~bo&CD=Ti&3W2W``>fIUZUZth2jZ-;-b8^n}N1*ZHuIzmH3wZuLW z!_S+~&1;<((o%*0{IR(5L1zVW=Shu>xwLC}(XAe zg1JBZY}O0k-)UiaP!Tcz$B!Q+o#%sO>Kq@Bc_|ybD|VPF%Yg1T?cA=tUbC zcbbt6b^wu}L<78_0N&cfX5DF~cg4&rqTf{`^IKlDj)X5gV7pJUKVM(fN}iM_AuB2< z%HbV~dT4L+9`L!rz4c|vXOEeSjif5B;J$fIiu4upXS1B4zwWl+gWL6ph4O}(vT3qm z{NJv(RwKEt^KhkbaqK2jzm2#Wx%-^%plEL})khTN-W2}SGG$DClXnRDO!E!E0d);B zhdj67r-Q3MrGAu|(_^UkWSC?FC$B6NFUQQh7ed$Qe79%wTLp0rZG=Bebl1MLV0-_H z|L)}C!_QL%c-k$9C~AGxJ_cff=F7pBbHyJw2dYEFSw~`xd=@_qhA+b@EU@gcaWfc(QDN*tt$uK$;--C zk{lcec%6O*I$xWXb);`1%`#VPdjEn8Y{oHt7gq4QMnKl+TM!BGT2Dgy9b;&9Nc7k8 zEYXQx*EA}~%9qLhHVxoRJo}_i*@I&TlFhgjYJ&3I4p%Uz9H^*uZ#kHF-xCzc`j277 z|EXX9qiOnoQBY!{m*e-Z>)_QzWAbm?S{5TpO2!WoMJPXcS$-Q^nAGKEiZh@18Sx7Iua#=HV3Q6mU3gi$nKE zSvL_xJ*qD!xP8k5FHO0tO}0DL#>@BP^~j8{#CxVGU225CWD@R4x&7ry~9|=P_u2)GlpkI&zNP~xv4Pem9zK^_5syq{3gFKAU}V5dPWpt)sw1P z{rK2iuAE(f&FU}W+L>*ezx3q_jSL*zM3G(gC!sv6QwZ^0+irLUGgFrp9&-C`S3$go z0=My-bhRS#!Te>II@c?NxR(!kjp}d|2Ve9M#ZqLq1XOu_G_ZCilJ66lF&Eqi{tCx+^y52 zb{{TW_VNnQB5Y6CUV!ew4vx+wVXyMe&|9bkWAeIEH{e;$hJ6zyTM_HBSa$BIiLs`h zQ|iE1jq)EF+8`gZp2|327;5hGn8YcLNa!8r--b^;Z#h;` zR#t|?n9XRj6WW~kBg@unS)w|AV;n_M=;Hk4rAD6dv6%|5)f`TjlgqP%Eaa%so ziHDoqZL_m+s}e(F(*}9=2**d&tYo{G-)SKewjH@0>o0j2XG0dkkZEP2$t8+|^w||0@bl#H7*%VC9 ztbmpplCumycRiG^%gsC8jtwjdK|0^}PXDm`@8llM182c3-~S-^5Ke->v3jrkhP0O8 zR#g)`9{JY4)!(XDS0m3sFTP}}_@jLKDuArNLD@ugtZc#R0GCn2ZvB@lqtq>7JJ>GN zoV0_9#S;}ED??umTall-pZo^rvv0$i$y?rWiR=>H3e@7*vt&83dUpF#1jCc+`kc*X z<>TtQ6IGBT-up@*g^Vs>;3?wmqoM2G(M0wTE$!XR508v17duG*PU`tLQUY|ID2wCO zkKID*CVGo0ieK`|TngETl!FO;Wi&^_Co)=A`pG{yy+V!|peL6PV7oLxuOLUmz7=FK zvcMUAn9iXw6YwW5R*Xx)B~o6B_{oTQ#p0+?hz&xx(5#We50a3uh#*rCVrf~6^>6!D zH@Yt|1AXkU^RIgD)PcGo-3C^QS+rtta$t&#>woc(Ng)J3RgmD%-;`W%=Fh;kibXkK zXo;{yrL_CxI+pvrscyMxtHmpJI5qV=VWv}WhQ|4a%G!!wKDhg>@Oyop+__4H_+y7j z&uvDHCMH^K#pp?7z#Wqvp`^EoZ2>C<50~T(a5Y(GJJpVQ&ddhn>L=xf}yRhk{9oqHt15Mq4!zTHFq1+ zOs`|@v#3V8vQ+XZn;GNNyz#Y2EW!Wh&dYz%(*NIpUYH)+?#jxFO8mvNWwrX77lMe0 z+2nee`dHWe>nPsJzxj|RJovX88yoO{W7AmXk5!seU9HQma)8={;9~jDxD8^Y7!RzU zXl=dH78P|Qc%}4-ySMB>Ui50eu@Sb{KQQq6G^v*X>-kB^doCG7n_h$K|9Dcdz7lDYUCPD!?#q$+yyB2a+_%OkU^O!p#bLt&Qvig=6B(>}( z>ZJ{mOw#(AWfKXXQd{4jdJGywiT1!9MahuO*y2pzyqE9fRC}h@<^*$GXaEF2WU4aa zFr&EuI|};;JIJakz{C%vwp4ITDQMVf=CjC3`bq!mDrX&&QIN81` z#^vS#7tpiwWJaFbKLfJf`=<4(__8n3lamV?tZR@fuDoQ^wSj>f(Jd}NZMC+(wkj}p z&qrou-S#~<6`-96B4fQyp#CYsQPp{>p=K1hA!2an3q43*voLLnx-xhuUOc5piTPDQ zm5!Z~^VP>xLd%z4#Mx&$h8mZzg-s!JAv)wXxTIBg0QBTgtFppq>kzLB)q(_t80MT; z;V<_X&-CU$);HgrLU+0=K!yUQ2Y|>xz0g;IHYfJyt-9?mryptUSltHwcLu04XK*G* zIg^x#CrHql7kxdCz8x?`PYR8CTORaf3O9ZC!J3h4!A<3A`sqhY()4%zvjM{1pwl%#QSJ2lL>dtc17548wSMpAq z)W+Tmn++3^!Ty*pO_|M-)j1F0=YR#_oD$o z)g`y@N}ss2Fse#z!~62#?*JMlMH#3S?<7Gky3zvsJ+Aha;m}on4VN|Po~~D5GSXU; z8DB&1yOpgdc@wK4YUY4Cv{1KyfA7DTf~HRV)W*%wVmWB2a9*zXVEc}tBkC5?$@SRI z*SAYdjP4fECvn-W(0R^P{=yJBkfw=JW)yUaxROENs2(4EL?GVeIcdtgmWMZKC3aT( z^%a{n4v=QK%zt4%Cijh<)%+y)nAc4isUxuf^C;0W1$;m-8F0v$xT!52s+#c3ET)>g zvkJ)I>iWn<4;t$1U6?!YLNz?N&@SF`F^viYYS6LEIbblqcgDTA&}@S=EnQBR_JE=< z;sQ}RiG$iIoNT*WXLHwF`)OtOPKQ41uqTqJ(3`s&zR?Aup+C=zLwj~oG||;n{cNKV z>Fd4u2`oK(4niH;nVew$0JH}pPrU6&v?S2JhWSmNJnux8C&4uk5S_2G)Dz3UdoifW z+Z8%GezQ}|eHw{;NRICIiNHc@)k~Z7FKX5S%o)pwVDNZQ?h1Id-lgI%$}AkOIJp8j z&0U(svegk?UQ*Gb$I1$8(o+>>BU@m@OUIL+L6;zQ!7EoClJOvF@&__xR3&ppuqC+Y zB4c+`C|6kEA^+T=9h|wNN(@hEK5?*OHtfcXFslbc39C+dRqKdy{w}_lun#*i|8QsG zYP_bYH%u^-rNlI=hj+UZ;BMnf6tceFe&|`cy!e~_ zMf^t|xQb@cja8^Gv01%TRlj0;n|&~WmcNOUQ(!_DgIv|HWA)L!eB6DH=EDPs*CfO2^?-zN=(1wgrAJ_} z(oeqd=a-5<;FlS+BqSt*Ek2bsEwf*5Mx!r3egZ9>Bd~pV#^nmrw7iUUv2OYX$XI7` zO4g=d(q$5YI)s07FgQH#?%V6>?tj+H@b9t>PLAEfqUIVJhn??>WN&Ysn%+jae2;&a z_#P3g;}Ub?hv;+Y+d19slH@)=jLeS{yItRX`%{YVg^=z47N|Q4-Td7qGJCU%(QNA9~xLUnFcZ(SLze;4T3BcCfGUNB+m{0J%UNei;Ghe6Fe#_&1{e^ zz;f$>n$k`96;N3j_x5WmVLSu``26KRVwYmdjMV2|f7I=WQHD1LZQb=xA&e(KegKNDWo=MM@5AdB%V zLRKVuu`;c{wF3I$^CJ-8*+T&Io;{l<3?`9gm)%}bb`RgLLq3ROg*>>Zj%xB|l|Zvi zJ}6^V>HK>WOYdbO9E2MaEO4&-jm;c)osUT+!i>$zK+`x450G07yz|HC+cmzKvY|$r z&7Bk;Nc>CNT!pC!BoTt)P$gVvr7lOl+@Ba}o^JtE=hY7(Uw}B$fbN2*<7)2zPh+$D dq^q~Ma7(2q?rE18Y$S}Mt*)