Use icon instead of self-built triangle array

triangle array was redundant

src/refine_surface.f90

@@ -77,7 +77,6 @@ type (edge), dimension(:), allocatable :: edswap
 integer,dimension(:)  ,allocatable :: nedgepernode
 integer,dimension(:)  ,allocatable :: refine_list,remove_list
 integer,dimension(:,:),allocatable :: nodenodenumber,nodeedgenumber
-integer,dimension(:,:),allocatable :: triangle
 logical,dimension(:),allocatable::refine,remove,removep,blocked
 double precision,dimension(:),allocatable::edlen,edlenp
 character*72  :: shift
@@ -86,7 +85,7 @@ integer tr1,tr2
 integer ntriangle
 integer nedge
 integer nedgen,nsurfacen
-double precision dist0,distmax,prod,distmin,eps
+double precision dist0,distmax,distmin,eps
 double precision maxdist1,maxdist2
 double precision, external :: dist
 
@@ -140,6 +139,8 @@ do ie=1,ntriangle !loop over all triangles
     inodepp=surface%icon(kpp,ie)
 
     !get edge information of right hand triangle
+
+    !is this correct? nedgepernode is 0 in the first run, so the loop will execute with j=1 and j=0 at the start...
     do j=1,nedgepernode(inodep)
       if (nodenodenumber(j,inodep).eq.inode) then
         jedge=nodeedgenumber(j,inodep)
@@ -172,44 +173,6 @@ enddo
 
 if (iedge/=nedge) call stop_run ('pb_b in refine_surface$')
 
-!Create index of edge numbers for all triangles
-!Edge numbers are sorted by number
-
-allocate (triangle(ntriangle,3),stat=err)
-triangle=0
-
-do iedge=1,nedge
-    !left-hand triangle
-    tr1 = ed(iedge)%t1
-    if (tr1 /= 0) then
-        !fill edge information incremetally
-        if (triangle(tr1,1)==0) then
-            triangle(tr1,1)=iedge
-        elseif (triangle(tr1,2)==0) then
-            triangle(tr1,2)=iedge
-        elseif (triangle(tr1,3)==0) then
-            triangle(tr1,3)=iedge
-        else
-            !problem if we get here
-            call stop_run ('Too many edges to one triangle')
-        end if
-    end if
-    !right-hand triangle
-    tr2 = ed(iedge)%t2
-    if (tr2 /= 0) then
-        !fill edge information incremetally
-        if (triangle(tr2,1)==0) then
-            triangle(tr2,1)=iedge
-        elseif (triangle(tr2,2)==0) then
-            triangle(tr2,2)=iedge
-        elseif (triangle(tr2,3)==0) then
-            triangle(tr2,3)=iedge
-        else
-            !problem if we get here
-            call stop_run ('Too many edges to one triangle')
-        end if
-    end if
-end do
 
 !------------------------------------
 ! Calculate length for all edges
@@ -257,20 +220,24 @@ do iedge=1,nedge
        dist0 = edlen(iedge)
        !Is this edge longer than distmax?
        if (dist0>distmax) then
+           !maxdist in right hand triangle
            tr1 = ed(iedge)%t1
-           maxdist1 = max( edlen(triangle(tr1,1)),edlen(triangle(tr1,2)),edlen(triangle(tr1,3)) )-eps
+           maxdist1 = max( edlen(surface%icon(1,tr1)), & 
+                           edlen(surface%icon(2,tr1)), &
+                           edlen(surface%icon(3,tr1)) )-eps
+           !maxdist in left hand triangle
            tr2 = ed(iedge)%t2
-           maxdist2 = max( edlen(triangle(tr2,1)),edlen(triangle(tr2,2)),edlen(triangle(tr2,3)) )-eps
+           maxdist2 = max( edlen(surface%icon(1,tr2)), &
+                           edlen(surface%icon(2,tr2)), &
+                           edlen(surface%icon(3,tr2)) )-eps
            !Is this edge one of the longest?
            if (dist0>=maxdist1 .OR. dist0>=maxdist2) then
                refine(iedge) = .true.
-               !block other edges in triagles, actual edge number does not matter:
-               !n<=iedge will not be visited in this loop again.
-               !
-               blocked(triangle(tr1,2))=.true.
-               blocked(triangle(tr1,3))=.true.
-               blocked(triangle(tr2,2))=.true.
-               blocked(triangle(tr2,3))=.true.
+               !block all edges in both triangles
+               do k = 1,3
+                   blocked(surface%icon(k,tr1))=.true.
+                   blocked(surface%icon(k,tr2))=.true.
+               end do
            end if
        end if
    end if
@@ -292,7 +259,6 @@ endif
 
 
 deallocate (edlen)
-deallocate (triangle)
 
 
 !-----------------------
@@ -338,12 +304,12 @@ if (params%remove_surf_pts) then
       enddo
     enddo
   endif
-  if (irem/=nrem) call stop_run ('pb_b in refine_surface$')
+  if (irem/=nrem) call stop_run ('pb_b in refine_surface')
   deallocate (remove)
   deallocate (removep)
 endif
 
-if (iadd/=nadd) call stop_run ('pb_b in refine_surface$')
+if (iadd/=nadd) call stop_run ('pb_b in refine_surface')
 
 deallocate (refine)
 deallocate (blocked)