Added quadratic shape function changes for computing node distances -

raulrn · Jul 15, 2013 · e4e504b · e4e504b
1 parent cc57541
commit e4e504b
Show file tree

Hide file tree

Showing 3 changed files with 43 additions and 16 deletions.
diff --git a/mpm/src/mpmutils_c.pyx b/mpm/src/mpmutils_c.pyx
@@ -91,6 +91,20 @@ def divergence( np.ndarray[ITYPE_t, ndim=2] cIdx,
             gg[ixc,1] -= pp[ii,1,0] * cg0 + pp[ii,1,1] * cg1; 
     return 0      
 
+def gradscalar( np.ndarray[ITYPE_t, ndim=2] cIdx, 
+              np.ndarray[FTYPE_t, ndim=3] cGrad, 
+              np.ndarray[FTYPE_t, ndim=2] pp,
+              np.ndarray[FTYPE_t, ndim=2] gg ):
+    # Send gradient of scalar particle field to grid
+    cdef int ii, jj, ixc
+    cdef int nParts = pp.shape[0]
+    cdef int nContrib = cIdx.shape[1]
+    for ii in range(nParts):
+        for jj in range(nContrib):
+            ixc = cIdx[ii,jj];  
+            gg[ixc,0] += pp[ii,0] * cGrad[ii,jj,0];
+            gg[ixc,1] += pp[ii,0] * cGrad[ii,jj,1];
+    return 0      
 
 def dotAdd( np.ndarray[FTYPE_t, ndim=3] pp,
             np.ndarray[FTYPE_t, ndim=3] qq ):

diff --git a/mpm/src/shapes/quad2.py b/mpm/src/shapes/quad2.py
@@ -30,23 +30,28 @@ def updateContribList( dw, patch, dwi ):
     # Update node contribution list
     nx = patch.Nc[0]
     h = patch.dX
+    hm = min(h)
     idxs = [0,1,2,nx,nx+1,nx+2,2*nx,2*nx+1,2*nx+2]
     S = np.zeros(h.size)
     G = np.zeros(h.size)	
 
     cIdx,cW,cGrad = dw.getMult( ['cIdx','cW','cGrad'], dwi )
     px,gx = dw.getMult( ['px','gx'], dwi )
+    gDist = dw.get( 'gDist', dwi )
 
     for ii in range(len(pVol)):
         cc = getCell( patch, px[ii] )	           
 
         for jj in range(9):	
             idx = idxs[jj] + cc 
-            r = px[ii] - gx[idx]	
+            r = px[ii] - gx[idx]
+            d = np.linalg.norm(r)
 
             for kk in range(len(r)):
                 S[kk],G[kk] = uSG( r[kk], h[kk] )
 
             cIdx[ii][jj] = idx
             cW[ii][jj] = S[0]*S[1]
-            cGrad[ii][jj] = G * S[::-1] 
+            cGrad[ii][jj] = G * S[::-1] 
+            gDist[idx] = max(0,gDist[idx], (1. - d/hm) )
+
diff --git a/mpm/src/shapes/quad2_c.pyx b/mpm/src/shapes/quad2_c.pyx
@@ -26,9 +26,9 @@ def updateContribList( dw, patch, dwi ):
     ng = patch.nGhost
     inpf = np.array([h,patch.X0, patch.dX])
     idxs = np.array([0,1,2,nx,nx+1,nx+2,2*nx,2*nx+1,2*nx+2])
-    labels = ['px','gx','cIdx','cW','cGrad']
-    px,gx,cIdx,cW,cGrad = dw.getMult(labels,dwi)
-    updateContribs( inpf, idxs, ng, th, px, gx, cIdx, cW, cGrad )
+    labels = ['px','gx','cIdx','cW','cGrad','gDist']
+    px,gx,cIdx,cW,cGrad, gDist = dw.getMult(labels,dwi)
+    updateContribs( inpf, idxs, ng, th, px, gx, cIdx, cW, cGrad, gDist )
 
 
 #===============================================================================        
@@ -38,45 +38,52 @@ def updateContribs( np.ndarray[FTYPE_t, ndim=2] inpf,
                     np.ndarray[FTYPE_t, ndim=2] gx, 
                     np.ndarray[ITYPE_t, ndim=2] cIdx, 
                     np.ndarray[FTYPE_t, ndim=2] cW, 
-                    np.ndarray[FTYPE_t, ndim=3] cGrad ):
+                    np.ndarray[FTYPE_t, ndim=3] cGrad,
+		    np.ndarray[FTYPE_t, ndim=1] gDist):
     # inpf - float input vector - h, dxdy, patch.X0, patch.dX
     cdef int nParts = px.shape[0]
     cdef int cc, idx 
     cdef int ii, jj, kk
-    cdef double x, r, h, sgn
+    cdef double r, h, sgn, hm
     cdef double* hh = [inpf[0,0],inpf[0,1]]
+    cdef double* x = [0.,0.]
     cdef double* pp = [0.,0.]
     cdef double* S = [0.,0.]
     cdef double* G = [0.,0.]
     cdef int* cix = [0,0]
-
+
+    hm = min(hh[0],hh[1])    
+    hm = sqrt(hh[0]*hh[0]+hh[1]*hh[1])
+
     for ii in range(nParts):
 
 	# Get Cell
         for kk in range(2):
             pp[kk] = px[ii,kk]
 
-            x = (pp[kk] - inpf[1,kk])/inpf[2,kk] + ng;
-            cix[kk] = int(floor(x))
-            if (x - 1.*cix[kk]) < 0.5:    cix[kk] -= 1
+            x[kk] = (pp[kk] - inpf[1,kk])/inpf[2,kk] + ng;
+            cix[kk] = int(floor(x[kk]))
+            if (x[kk] - 1.*cix[kk]) < 0.5:    cix[kk] -= 1
 
         cc = int(cix[1]*idxs[3] + cix[0]);
 
         for jj in range(9):
             idx = cc + idxs[jj];
+            x[0] = pp[0]-gx[idx,0]
+            x[1] = pp[1]-gx[idx,1]
+            d = sqrt( x[0]*x[0] + x[1]*x[1] )
 
             for kk in range(2):
-                x = pp[kk] - gx[idx,kk]
-                r = fabs(x)
+                r = fabs(x[kk])
                 h = hh[kk]
-                sgn = copysign(1.,x)
+                sgn = copysign(1.,x[kk])
 
                 if ( r < 0.5*h ):
                     S[kk] = -r*r/(h*h) + 3./4
-                    G[kk] = -2.*x/(h*h)
+                    G[kk] = -2.*x[kk]/(h*h)
                 elif ( r < 1.5*h ): 
                     S[kk] = r*r/(2.*h*h) - 3.*r/(2.*h) + 9./8
-                    G[kk] = x/(h*h) - sgn*3./(2*h)        
+                    G[kk] = x[kk]/(h*h) - sgn*3./(2*h)        
                 else: 
                     S[kk] = 0.
                     G[kk] = 0.	
@@ -85,5 +92,6 @@ def updateContribs( np.ndarray[FTYPE_t, ndim=2] inpf,
             cW[ii,jj] = S[0]*S[1]
             cGrad[ii,jj,0] = S[1]*G[0]
             cGrad[ii,jj,1] = S[0]*G[1]
+            gDist[idx] = max(0, max(gDist[idx], (1. - d/hm)))
 
     return 0