From 12dad7c28ecb27ea08a93f5fac9cefa57e48d218 Mon Sep 17 00:00:00 2001 From: sbosse <sbosse@uni-bremen.de> Date: Wed, 28 Aug 2024 21:43:15 +0200 Subject: [PATCH] Wed 28 Aug 21:38:52 CEST 2024 --- src/SimNDT/engine/EFIT2Dcython.pyx | 141 +++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) create mode 100644 src/SimNDT/engine/EFIT2Dcython.pyx diff --git a/src/SimNDT/engine/EFIT2Dcython.pyx b/src/SimNDT/engine/EFIT2Dcython.pyx new file mode 100644 index 0000000..3451a8d --- /dev/null +++ b/src/SimNDT/engine/EFIT2Dcython.pyx @@ -0,0 +1,141 @@ +cimport cython +cimport numpy as np + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +def sourceVel(np.float32_t[:,:] Vx, np.float32_t[:,:] Vy, + np.int32_t[:] XL, + np.int32_t[:] YL, + double source, int Stype, win, double dtdxx): + + cdef unsigned int i + cdef int ii, jj + cdef int N = XL.shape[0] + + if Stype == 0: + for i in range(N): + ii = XL[i] + jj = YL[i] + Vx[ii, jj] -= source*dtdxx * win[i] + + elif Stype ==1: + for i in range(N): + ii = XL[i] + jj = YL[i] + Vy[ii,jj] -= source*dtdxx * win[i] + + elif Stype ==2: + for i in range(N): + ii = XL[i] + jj = YL[i] + Vx[ii,jj] -= source* dtdxx * win[i] + Vy[ii,jj] -= source* dtdxx * win[i] + + + return (Vx, Vy) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +def sourceStress(np.float32_t[:,:] Txx, np.float32_t[:,:] Tyy, + np.int32_t[:] XL, np.int32_t[:] YL, + double source, int Stype, win, double dtdxx): + + cdef int N = XL.shape[0] + cdef unsigned int i + cdef int ii, jj + + if Stype ==0: + for i in range(0, N): + Txx[XL[i],YL[i]] -= source*dtdxx * win[i] + + elif Stype ==1: + for i in range(0, N): + Tyy[XL[i],YL[i]] -= source*dtdxx * win[i] + + elif Stype ==2: + for i in range(0, N): + Txx[XL[i],YL[i]] -= source* dtdxx * win[i] + Tyy[XL[i],YL[i]] -= source* dtdxx * win[i] + + + return (Txx, Tyy) + + +@cython.boundscheck(False) +@cython.wraparound(False) +def velocityVoigt(np.float32_t[:,:] Txx, np.float32_t[:,:] Txy, np.float32_t[:,:] Tyy, + np.float32_t[:,:] vx, np.float32_t[:,:] vy, np.float32_t[:,:] d_vx, np.float32_t[:,:] d_vy, + np.float32_t[:,:] BX, np.float32_t[:,:] BY, np.float32_t[:,:] ABS, double ddx, double dt): + + + cdef int M = Txx.shape[0] + cdef int N = Txx.shape[1] + cdef int i + cdef int j + + + for i in range(0,M-1): + for j in range(1,N): + d_vx[i,j] = (BX[i,j] * ddx) * ( Txx[i+1,j] - Txx[i,j] + Txy[i,j] - Txy[i,j-1] ) + vx[i,j] += dt * d_vx[i,j] + + + for i in range(1,M): + for j in range(0, N-1): + d_vy[i,j] = (BY[i,j] * ddx) * ( Txy[i,j] - Txy[i-1,j] + Tyy[i,j+1] - Tyy[i,j] ) + vy[i,j] += dt * d_vy[i,j] + + for i in range(0,M): + for j in range(0,N): + vx[i,j] *= ABS[i,j] + vy[i,j] *= ABS[i,j] + + + + return (vx, vy, d_vx, d_vy) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +def stressVoigt(np.float32_t[:,:] Txx, np.float32_t[:,:] Txy, np.float32_t[:,:] Tyy, + np.float32_t[:,:] vx, np.float32_t[:,:] vy, np.float32_t[:,:] d_vx, np.float32_t[:,:] d_vy, + np.float32_t[:,:] C11, np.float32_t[:,:] C12, np.float32_t[:,:] C44, + np.float32_t[:,:] ETA_VS, np.float32_t[:,:] ETA_S, np.float32_t[:,:] ETA_SS, + np.float32_t[:,:] ABS, double dtx): + + + cdef int M = Txx.shape[0] + cdef int N = Txx.shape[1] + cdef int i + cdef int j + + for i in range(1,M): + for j in range(1,N): + Txx[i,j] += ( ( C11[i,j] * dtx )*(vx[i,j] - vx[i-1,j]) + + ( C12[i,j] * dtx )*(vy[i,j] - vy[i,j-1]) + + ( (ETA_VS[i,j]) * dtx) * (d_vx[i,j] - d_vx[i-1,j]) + + ( (ETA_S [i,j]) * dtx) * (d_vy[i,j] - d_vy[i,j-1]) ) + + Tyy[i,j] += ( ( C12[i,j]* dtx )*(vx[i,j] - vx[i-1,j]) + + ( C11[i,j]* dtx )*(vy[i,j] - vy[i,j-1]) + + ( (ETA_VS[i,j])*dtx) *(d_vy[i,j] - d_vy[i,j-1]) + + ( (ETA_S [i,j])*dtx) * (d_vx[i,j] - d_vx[i-1,j]) ) + + for i in range(0, M-1): + for j in range(0, N-1): + Txy[i,j] += ( ( C44[i,j] * dtx )*(vx[i,j+1] - vx[i,j] + vy[i+1,j] - vy[i,j]) + + ( ETA_SS[i,j] * dtx )*(d_vx[i,j+1] - d_vx[i,j] + d_vy[i+1,j] - d_vy[i,j] ) ) + + for i in range(0,M): + for j in range(0,N): + Txx[i,j] *= ABS[i,j]; + Tyy[i,j] *= ABS[i,j]; + Txy[i,j] *= ABS[i,j]; + + + return (Txx, Tyy, Txy) + + + \ No newline at end of file