rm(list=ls())
#
#-------------------------------------------------------------------------------
# In diesem Skript werden die Komponenten des Spannungstensors dargestellt.
# Zu einer Ebene in Dreipunktedarstellung wird der Normaleneinheitsvektor
# berechnet und eingezeichnet. Aus dem Spannungstensor und dem Normaleneinheits-
# vektor wird der Spannungsvektor in dieser Ebene bestimmt und in die 
# Normal- und Scherspannung zerlegt.
# Aufgabe: Variieren Sie die Lage der Ebene und versuchen Sie eine Hauptrichtung
#          zu finden.
#-------------------------------------------------------------------------------
#
# Das ist der Spannungstensor:
S <- 0.4*matrix(
            c(-4.0, 3.0 , 2.0,
              3.0, 2.0 , 7.0,
              2.0, 7.0 , 2.2),nrow=3,byrow=T)
# Test auf Symmetrie
if(S[2,1]!=S[1,2]||S[3,1]!=S[1,3]||S[3,2]!=S[2,3])
stop("\n*** tensor not symmetrical ***")
print(S)
#-------------------------------------------------------------------------------
# Punkte zur Festlegung der Ebene im Raum (Dreieck):
#-------------------------------------------------------------------------------
p1 <- c(5,0,0)
p2 <- c(0,0,4)
p3 <- c(0,5,0)
#-------------------------------------------------------------------------------
#    Der Normaleneinheitsvektor der Ebene:
#-------------------------------------------------------------------------------
dp12 <- p2 - p1
dp13 <- p3 - p1
#
# Kreuzprodukt und Normalisierung:
nev <- c(dp12[2]*dp13[3] - dp12[3]*dp13[2],
         dp12[3]*dp13[1] - dp12[1]*dp13[3],
         dp12[1]*dp13[2] - dp12[2]*dp13[1])
nev0 <- -nev/sqrt(nev%*%nev) # Normaleneinheitsvektor
print(nev0)         
#-------------------------------------------------------------------------------
# Berechnung des Spannungsvektors
#-------------------------------------------------------------------------------
sv <- t(S%*%nev0)
print(sv)
#
# Berechung des Dreiecksschwerpunktes zum Einzeichen des Spannungsvektors:
#
x0 <- (1/3)*(c(p1[1]+p2[1]+p3[1],p1[2]+p2[2]+p3[2],p1[3]+p2[3]+p3[3]))
x1 <- x0 + sv # Endpunkt des Spannungsvektors
#
# Normalspannung:
sigma <- (sv%*%nev0)*nev0
#
# Scherspannung:
tau <- sv - sigma
#-------------------------------------------------------------------------------
# ab jetzt: "nur" noch Zeichnen der Ebene und der Pfeile
#-------------------------------------------------------------------------------
#
u <- c(0,0.001); v <- c(0,0.001) 
# Ein winziges Rechteck muss leider gezeichnet werden.
# Ich hoffe man sieht es nicht.
datmat <- outer(u,v)
#
persp(u,v,datmat, xlim = c(0,6), ylim= c(0,6), zlim = c(0,6),
     theta = 30, phi = 20, expand = 0.7, col = rgb(0,1,1),
     ltheta = 120, shade = 0.75, ticktype = "detailed",
     xlab = "X-Value", ylab = "\nY-Value", zlab = "Z-Value", 
     main="Spannungstensor und\nSpannungsvektor",r=4,border=NA) -> rotmat
#
# Einzeichnen der Dreiecke in den Koordinaten-Ebenen
#
x <- c(p1[1],p2[1],p3[1])
y <- c(p1[2],p2[2],p3[2])
z <- c(    0,    0,    0)
polygon (trans3d(x, y, z, pmat = rotmat), lwd=1, col = "gray90", border = "darkgreen",lty=1)
# 
# yz-Ebene
x <- c(0    ,    0,    0)
y <- c(p1[2],p2[2],p3[2])
z <- c(p1[3],p2[3],p3[3])
polygon (trans3d(x, y, z, pmat = rotmat), lwd=1, col = "gray70", border = "darkgreen",lty=1)
#
x <- c(p1[1],p2[1],p3[1])
y <- c(p1[2],p2[2],p3[2])
z <- c(p1[3],p2[3],p3[3])
polygon (trans3d(x, y, z, pmat = rotmat), lwd=2, col = "steelblue2", border = "darkgreen",lty=1)
#
# und hier die Pfeile       
source("script_bm_fun_Pfeile.r")
ep.tau    <- as.vector(x0 + tau)
ep.sigma  <- as.vector(x0 + sigma)
ep.sv     <- as.vector(x1)
ep.nev0   <- as.vector(x0+nev0)
# cat("\n*** hier sammer ***\n")
pfeil(x0,ep.sv,0.3,0.1,rotmat,"black")    
pfeil(x0,ep.sigma,0.3,0.1,rotmat,"blue")
pfeil(x0,ep.tau,0.3,0.1,rotmat,"blue")
pfeil(x0,ep.nev0,0.3,0.1,rotmat,"green")
text (trans3d(x1[1]-0.2, x1[2]+0.2, x1[3]+0.2, pmat = rotmat), "sv")
x <- c(ep.tau[1],ep.sv[1],ep.sigma[1])
y <- c(ep.tau[2],ep.sv[2],ep.sigma[2])
z <- c(ep.tau[3],ep.sv[3],ep.sigma[3])
lines(trans3d(x,y,z, pmat=rotmat), lty=2)
#-------------------------------------------------------------------------------
# und jetzt: die Komponenten des Spannungstensors
#-------------------------------------------------------------------------------
#
# die Startpunkte der Pfeile:
#
s0.xy <- (p1+p3)/3
s0.xz <- (p1+p2)/3
s0.yz <- (p2+p3)/3
#
# die Pfeilchen in der yz-Ebene:
# 
sigma.x <- c(S[1,1],     0,     0)
tau.xy  <- c(0     ,S[2,1],     0)
tau.xz  <- c(0     ,     0,S[3,1])
pfeil(s0.yz,s0.yz+sigma.x,0.3,0.1,rotmat,"gray80")
pfeil(s0.yz,s0.yz+tau.xy,0.3,0.1,rotmat,"gray80")
pfeil(s0.yz,s0.yz+tau.xz,0.3,0.1,rotmat,"gray80")
#
# die Pfeilchen in der xz-Ebene:
# 
sigma.y <- c(     0,S[2,2],     0)
tau.yx  <- c(S[1,2],     0,     0)
tau.yz  <- c(0     ,     0,S[3,2])
pfeil(s0.xz,s0.xz-sigma.y,0.3,0.1,rotmat,"black")
pfeil(s0.xz,s0.xz-tau.yx,0.3,0.1,rotmat,"black")
pfeil(s0.xz,s0.xz-tau.yz,0.3,0.1,rotmat,"black")
#
# die Pfeilchen in der xy-Ebene:
# 
sigma.z <- c(     0,     0,S[3,3])
tau.zx  <- c(S[3,1],     0,     0)
tau.zy  <- c(0     ,S[3,2],     0)
pfeil(s0.xy,s0.xy-sigma.z,0.3,0.1,rotmat,"gray")
pfeil(s0.xy,s0.xy-tau.zx,0.3,0.1,rotmat,"gray")
pfeil(s0.xy,s0.xy-tau.zy,0.3,0.1,rotmat,"gray")
#
# Einzeichnen der Startpunkte:
#
x <- c(s0.xy[1],s0.xz[1],s0.yz[1])
y <- c(s0.xy[2],s0.xz[2],s0.yz[2])
z <- c(s0.xy[3],s0.xz[3],s0.yz[3])
points(trans3d(x,y,z, pmat=rotmat), pch=20,col="blue")
#
#stop("\n*** uff ***") 
#
#   --- Schluss ---
#