diff --git a/LoopStructural/interpolators/cython/dsi_helper.pyx b/LoopStructural/interpolators/cython/dsi_helper.pyx index 274e60748..63a158586 100644 --- a/LoopStructural/interpolators/cython/dsi_helper.pyx +++ b/LoopStructural/interpolators/cython/dsi_helper.pyx @@ -25,6 +25,7 @@ def cg(double [:,:,:] EG, long long [:,:] neighbours, long long [:,:] elements,d cdef double [:,:] e1 cdef double [:,:] e2 cdef double area = 0 + cdef double length cdef long long [:] idl = np.zeros(4,dtype=np.int64) cdef long long [:] idr = np.zeros(4,dtype=np.int64) for e in range(ne): @@ -67,10 +68,13 @@ def cg(double [:,:,:] EG, long long [:,:] neighbours, long long [:,:] elements,d norm[1] = v1[2]*v2[0] - v1[0]*v2[2] norm[2] = v1[0]*v2[1] - v1[1]*v2[0] + length = np.linalg.norm(norm) # we want to weight the cg by the area of the shared face # area of triangle is half area of parallelogram # https://math.stackexchange.com/questions/128991/how-to-calculate-the-area-of-a-3d-triangle - area = 0.5*sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2])#np.linalg.norm(norm) + area = 0.5*length#sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2])#np.linalg.norm(norm) + for i in range(3): + norm[i]/=length for itr_left in range(Na): idc[ncons,itr_left] = idl[itr_left] for i in range(3):