From ff9b242076abbaa479069d2ee26296c1b8e281a6 Mon Sep 17 00:00:00 2001
From: Hoffmann <hoffmann@mathematik.tu-darmstadt.de>
Date: Fri, 8 Nov 2024 11:23:39 +0100
Subject: [PATCH] changes in FV scheme

---
 2d/myfun.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/2d/myfun.py b/2d/myfun.py
index 11ae8ac..d16221f 100644
--- a/2d/myfun.py
+++ b/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))
-- 
GitLab