changes in FV scheme

2d/myfun.py

@@ -1161,8 +1161,8 @@ def finitevolumescheme_rho(u_old,ht,numtri,tri,K,eps,v,f) :
                     [vKE,areaE] = getvKE(K[i],E,v) # avrg over v*nKE ( = 0 if E on boundary due to neumann boundary w.r.t. c )
                     dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of triangles
 
-                    if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case centered flux to avoid divide by 0 or log(0)
-                        sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i]) + vKE*areaE*1/2*(u[i]+u[neighbors[j]])
+                    if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case flux=0 to avoid divide by 0 or log(0)
+                        sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i])
                     else :
                         sum_EK[i] = sum_EK[i] - eps*areaE/dE*(u[neighbors[j]]-u[i]) + vKE*areaE*((u[i]-u[neighbors[j]])/(np.log(u[i])-np.log(u[neighbors[j]])))
 
@@ -1199,9 +1199,9 @@ def finitevolumescheme_rho(u_old,ht,numtri,tri,K,eps,v,f) :
                     [vKE,areaE] = getvKE(K[i],E,v) # avrg over v*nKE
                     dE = getdistE(K[i],K[neighbors[j]]) # get distance between circumcenters of triangles
 
-                    if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case centered flux to avoid divide by 0 or log(0)
-                        valK = valK + eps*(areaE/dE) + vKE*areaE*1/2
-                        valL = valL - eps*(areaE/dE) + vKE*areaE*1/2
+                    if np.abs(u[i]-u[neighbors[j]]) < 10**-6 or u[i] < 10**-6 or u[neighbors[j]] < 10**-6 : # in this case flux=0 to avoid divide by 0 or log(0)
+                        valK = valK + eps*(areaE/dE)
+                        valL = valL - eps*(areaE/dE)
                     else :
                         valK = valK + eps*(areaE/dE) + vKE*areaE*((np.log(u[i])-np.log(u[neighbors[j]])-1+u[neighbors[j]]/u[i])/((np.log(u[i])-np.log(u[neighbors[j]]))**2))
                         valL = valL - eps*(areaE/dE) + vKE*areaE*((-np.log(u[i])+np.log(u[neighbors[j]])-1+u[i]/u[neighbors[j]])/((np.log(u[i])-np.log(u[neighbors[j]]))**2))