First attempt at implementing the surface stabilization of Kaus et al., 2010

src/build_system_wsmp.f90

@@ -68,7 +68,7 @@ double precision weightel_glob(osolve%nleaves)
 
 integer, external :: ifind_loc_wsmp
 double precision,dimension(:,:),allocatable::lsf_el
-double precision ael(8*ndof,8*ndof),bel(8*ndof),patch1(5,5),patch2(3,3)
+double precision ael(8*ndof,8*ndof),bel(8*ndof),patch1(5,5),patch2(3,3),elem_dz
 double precision,dimension(:),allocatable::penalty_loc,penalty_glob,eviscosity,weightel
 double precision r0,s0,t0,rst,amin
 integer i,j,k,itodo,idof,i1,i2,iii,jjj,err,jj,i1loc,k1,k2,ii
@@ -172,6 +172,8 @@ do ileaves=1,osolve%nleaves
          ! transfer elemental matrix and rhs into wsmpmatrix structure
          ! different for symmetric or nonsymmetric matrices
 
+         do i=1,params%mpe
+
          do k1=1,params%mpe
             do idof1=1,ndof
                i1=(vo%ftr(osolve%icon(k1,ileaves))-1)*ndof+idof1
@@ -189,6 +191,17 @@ do ileaves=1,osolve%nleaves
                         endif
                      enddo
                   enddo
+                  if (kfix(ileaves) < 1) then
+                    if (idof1 == 3) then
+                      iloc=iffind_loc_wsmp(i1loc,i1,ja,ia,n_iproc,nz_loc)
+                      elem_dz=(osolve%icon(5,ileaves)-osolve%icon(1,ileaves))*params%vex
+                      if (k1 < 5) then                        
+                        avals(iloc)=avals(iloc)+(params%dt*weightel(ileaves))/(2*elem_dz)
+                      else
+                        avals(iloc)=avals(iloc)-(params%dt*weightel(ileaves))/(2*elem_dz)
+                      endif
+                    endif
+                  endif
                   b(i1loc,1)=b(i1loc,1)+bel(iii)
                endif
             enddo