DOUG 0.2

GeomInterp.f90

Go to the documentation of this file.
00001 ! DOUG - Domain decomposition On Unstructured Grids
00002 ! Copyright (C) 1998-2006 Faculty of Computer Science, University of Tartu and
00003 ! Department of Mathematics, University of Bath
00004 !
00005 ! This library is free software; you can redistribute it and/or
00006 ! modify it under the terms of the GNU Lesser General Public
00007 ! License as published by the Free Software Foundation; either
00008 ! version 2.1 of the License, or (at your option) any later version.
00009 !
00010 ! This library is distributed in the hope that it will be useful,
00011 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013 ! Lesser General Public License for more details.
00014 !
00015 ! You should have received a copy of the GNU Lesser General Public
00016 ! License along with this library; if not, write to the Free Software
00017 ! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
00018 ! or contact the authors (University of Tartu, Faculty of Computer Science, Chair
00019 ! of Distributed Systems, Liivi 2, 50409 Tartu, Estonia, http://dougdevel.org,
00020 ! mailto:info(at)dougdevel.org)
00021 
00022 module GeomInterp
00023 ! This module is intended for bi/trilinear interpolation matrix generation
00024 
00025 ! General idea:
00026 ! As in CreateRestrict, a top down approach is taken. The tree is walked
00027 ! a coarse element at a time from coarsest to finest. 
00028 ! Bi/trilinear interpolation assume we know values in the corners of a 
00029 ! square/cube. In the grid elements we do, but in the refinements, we might
00030 ! not always have that information. In such a case the corner values usually
00031 ! need to be interpolated linearly from values connected to that point.
00032 ! Parsing the tree in a top down fasion has the advantage that we always
00033 ! know what points and with what coeficients to use when interpolating
00034 ! into the corners of the previous refinement. That info can rather easily
00035 ! be used to find the values in the corners of the new and finer element.
00036 ! This mechanism guarantees we know how to get the values at the corners of 
00037 ! the element once we get to it. Since we also know the center point of the
00038 ! element and want to use that for interpolation, we also subdivide the element
00039 ! we get into 4/8 others before actual interpolation. 
00040 ! Now we could just use these values and do the interpolation for each fine
00041 ! point based on the coarse coefficients. The trouble with that is that it
00042 ! might be quite slow. To speed up the process, some simple math is used, 
00043 ! that allows us to bring the single point interpolation costs down to
00044 ! 4 multiplications and 3 additions in bilinear and 12 mult. 7 add. for trilin.
00045 ! (might be overoptimizing but the inefficiency with which it was done in
00046 ! previous doug left a deep mark on my psyche - author).
00047 
00048 ! Details:
00049 ! I will not go into the finer points of mathematics - they should be checkable
00050 ! by anyone with a descent understanding of school algebra and the essence of
00051 ! bi/trilinear interpolation (which can be found explained in many places in 
00052 ! the internet). The actual interpolation algorithm stores the element bounds
00053 ! and coefficients of all previous levels for future use, which ensures no
00054 ! values are calculated more than once. Other than the math, there is 
00055 ! relatively little complicated in here (although one might want to reread
00056 ! the commentary of CreateCoarse). Math however is nearly impossible to explain
00057 ! by text alone, especially in the case of geometry.
00058 ! There is one result that might need mentioning - namely, that 
00059 ! the introduction of hanging nodes does not increase the upper bound of
00060 ! coarse elements used for interpolation for one fine one. That bound is tight
00061 ! and in both cases is 2**M%nsd+max(C%mlvl-1,0). The proof is not impossible:
00062 ! One just has to note that nearly always out of four corner points, 
00063 ! 3 have a coarse node used that no other one uses and that property can be
00064 ! proven to hold in descent assuming a trick is used "adding" a discarded 
00065 ! unique node to one lacking it with a coefficient zero should it be needed.
00066 ! Again - it is hard to explain with only words, as mathematics tends to be.
00067 
00068 ! NOTE: I have left debug output in this code commented out instead of deleted.
00069 ! The reason is that this code might need some rewriting (either towards
00070 ! readability or towards efficiency) and is sadly quite easy to break in the
00071 ! process. These debug snippets might help the process along.
00072 
00073 contains
00074 
00075     !! Calculates the node prolongation matrix framework, then call the filler
00076     !! Later rebuild the filled node prolong mat into a freedom restrict mat.
00077     subroutine CreateGeneralRestrict(C,M,R,genmat)
00078         use RealKind
00079         use CoarseGrid_class
00080         use SpMtx_class
00081         use SpMtx_util
00082         use globals, only:stream
00083         
00084         implicit none
00085 
00086         !! Coarse Grid whose structure to use
00087         type(CoarseGrid), intent(inout) :: C
00088         !! Fine Mesh for which to build
00089         type(Mesh), intent(in) :: M
00090         !! The restriction matrix to be created
00091         type(SpMtx), intent(out) :: R
00092   
00093         !! gendata is used to create aux data of size getsize 
00094         !!  for interpolation from points in pts.
00095         !! Output has to be so that all the elements concerning 
00096         !!  one fine node are consecutive
00097         interface
00098             subroutine genmat(A,C,M)
00099                 use Mesh_class
00100                 use SpMtx_class
00101                 use CoarseGrid_class
00102                 
00103                 type(SpMtx),intent(inout) :: A
00104                 type(Mesh), intent(in) :: M
00105                 type(CoarseGrid), intent(in) :: C
00106             end subroutine genmat
00107         end interface
00108         
00109         type(SpMtx) :: NP ! node prolongation matrix - restrict is transposed
00110        
00111         integer :: cnt
00112         integer :: i, j, k, pn, nd, f
00113         integer :: tnsd
00114         real(kind=xyzk) :: pts(M%nsd,2**M%nsd+C%mlvl)
00115         integer :: pinds(2**M%nsd+C%mlvl)
00116         real(kind=xyzk),pointer :: pt(:)
00117 
00118         ! For the reverse of freemap
00119         integer :: flist(C%nlfc)
00120         integer :: fbeg(C%nct+1)
00121         integer :: faux(C%nct)
00122         integer :: lbounds(M%lnnode), ubounds(M%lnnode)
00123 
00124 
00125         tnsd=2**M%nsd
00126 
00127         !****************************************************
00128         ! Create the Node Prolongation matrix structure
00129         !****************************************************
00130 
00131         ! Calculate nnz        
00132 
00133         cnt=0
00134         
00135         ! Increment by initial grid
00136         do i=1, C%elnum
00137             cnt=cnt+tnsd*C%els(i)%nfs
00138         enddo
00139 
00140         ! Increment by refinements
00141         do i=1, C%refnum
00142             cnt = cnt + &
00143                      C%refels(i)%level*(C%refels(i)%lend-C%refels(i)%lbeg+1)
00144         enddo
00145 
00146         ! Allocate the node prolong matrix
00147         NP=SpMtx_newInit(cnt,nrows=M%lnnode,ncols=C%nct)
00148 
00149         !****************************************************
00150         ! Create the prolongation data for NP
00151         !****************************************************
00152         call genmat(NP,C,M)
00153       
00154 !        call SpMtx_printRaw(NP)
00155 
00156         !****************************************************
00157         ! Create the reverse of freemap
00158         !****************************************************
00159 
00160         ! Find the counts of freedoms associated with nodes
00161         faux=0;
00162         do i=1,C%nlfc
00163             faux(C%cfreemap(i))=faux(C%cfreemap(i))+1
00164         enddo
00165 
00166         ! Create fbegs
00167         fbeg(1)=1;
00168         do i=1,C%nct
00169             fbeg(i+1)=faux(i)+fbeg(i)
00170         enddo
00171         
00172         ! Create flist
00173         faux=fbeg(1:C%nlfc)
00174         do i=1,C%nlfc
00175             nd=C%cfreemap(i)
00176             flist(faux(nd))=i
00177             faux(nd)=faux(nd)+1
00178         enddo
00179 
00180         !****************************************************
00181         ! Create R based on NP
00182         !****************************************************
00183 
00184         ! Calculate the upper bound for nnz
00185         cnt=0
00186         do i=1,NP%nnz
00187             cnt=cnt+fbeg(NP%indj(i)+1)-fbeg(NP%indj(i))
00188         enddo
00189 
00190         ! Calculate l and ubounds
00191         lbounds=0; ubounds=0
00192         lbounds(NP%indi(1))=1
00193         do i=2,NP%nnz
00194            if (NP%indi(i-1)/=NP%indi(i)) then
00195                if (ubounds(NP%indi(i-1))/=0 .or. &
00196                    lbounds(NP%indi(i))/=0) &
00197                      call DOUG_abort("Node restriction matrix not consolidated")
00198                ubounds(NP%indi(i-1))=i-1
00199                lbounds(NP%indi(i))=i
00200            endif
00201         enddo
00202         ubounds(NP%indi(NP%nnz))=NP%nnz
00203 
00204         ! Allocate the matrix
00205         R=SpMtx_newInit(cnt,nrows=C%nlfc,ncols=M%nlf)
00206         allocate(R%M_bound(M%nlf+1))
00207 
00208         ! Fill it
00209         f=1;         
00210         do i=1,M%nlf    
00211             nd=M%lfreemap(i)
00212             R%M_bound(i)=f
00213             if (nd/=0) then
00214             do j=lbounds(nd),ubounds(nd)
00215                 do k=fbeg(NP%indj(j)),fbeg(NP%indj(j)+1)-1
00216                 if (NP%val(j)>eps) then
00217                     R%val(f)=NP%val(j)
00218                     R%indi(f)=flist(k)
00219                     R%indj(f)=i
00220                     f=f+1
00221                 endif
00222                 enddo
00223             enddo
00224             endif
00225         enddo
00226         R%M_bound(M%nlf+1)=f
00227 
00228         ! Set the matrix to be row-sorted format
00229         R%arrange_type=D_SpMtx_ARRNG_COLS
00230 
00231         ! Destroy the node version
00232         call SpMtx_Destroy(NP)
00233         
00234         ! Give back the unneeded memory in P
00235         if (R%nnz>f-1) call SpMtx_resize(R,f-1)
00236 
00237         if (sctls%verbose>3) &
00238              write (stream,*) "Restriction matrix has ",f-1," elements"
00239 
00240         
00241    end subroutine CreateGeneralRestrict
00242 
00243    ! Calculate the next coefficients - refinement with center ct towards pt
00244    subroutine CalcNextCoefs(ct,ci,hnds,minv,maxv,coefs,cinds,csz,nsd,pt,cfout,ciout, ncsz)
00245        use RealKind
00246        use CoarseGrid_class
00247 
00248        real(kind=xyzk),intent(in) :: ct(:) ! coordinates of the current center
00249        integer, intent(in) :: ci ! index of the current center for ciout
00250        integer, pointer :: hnds(:) ! hanging nodes for the current center
00251        real(kind=xyzk),intent(in) :: minv(:),maxv(:)
00252        real(kind=xyzk),intent(in) :: coefs(:,:)
00253        integer, intent(in) :: cinds(:)
00254        integer, intent(in) :: csz ! number of coarse grid els in use in coefs
00255        integer, intent(in) :: nsd
00256        real(kind=xyzk),intent(in) :: pt(:)
00257        real(kind=xyzk),intent(out) :: cfout(:,:)
00258        integer, intent(out) :: ciout(:)
00259        integer, intent(out) :: ncsz
00260 
00261        integer :: i, k, ncnt
00262        integer :: ni(2*nsd+1), nc(2*nsd+1) ! new node indices and corners
00263        real(kind=xyzk) :: m(nsd) ! interpolation multipliers
00264 
00265        integer :: dir,dh(nsd),drev,rev(nsd)
00266        
00267        if (.not.associated(hnds)) then
00268            allocate(hnds(2*nsd)); hnds=0
00269        endif
00270 
00271        ! Set new count to 1 (if no hanging nodes are used, it is valid)
00272        ncnt=1 
00273 
00274        ! calc the interpolation multipliers
00275        m=(maxv-ct)/(maxv-minv) ! Coefficients for the min values in each dir
00276 
00277        ! Fill the outbound structures with initial values
00278        ncsz=csz; ciout(1:csz)=cinds(1:csz)
00279 
00280 !----------------------------------------------------------------------
00281 
00282        ! Find the direction of the element and the array to reverse it
00283        dir=1; drev=1; k=1
00284        do i=1,nsd
00285            if (pt(i)<ct(i)) then
00286                dir=dir+k; rev(i)=-k; dh(i)=nsd+i
00287            else
00288                drev=drev+k; rev(i)=k; dh(i)=i
00289            endif
00290            k=2*k
00291        enddo
00292        ! drev = <index of opposite corner to dir> = dir+sum(rev)
00293 
00294        ! The fixed point
00295        cfout(dir,1:csz)=coefs(dir,1:csz)
00296 
00297        ! The new center of refinement used
00298        cfout(drev,1:csz)=0.0_xyzk; nc(1)=drev; ni(1)=ci ! Add the new center
00299        
00300        if (nsd==2) then
00301            ! The nodes directly connected to the new center
00302            ! To do it for 3d case, one should use bilinear interpolation
00303            !  on the faces of the cube
00304            do i=1,nsd
00305                 if (hnds(dh(i))>0) then
00306                    cfout(drev-rev(i),1:csz)=0.0_xyzk; ncnt=ncnt+1
00307                    nc(ncnt)=drev-rev(i); ni(ncnt)=hnds(dh(i))
00308                 else
00309                    if (rev(i)>0) then ! Fixed node is max
00310                        cfout(drev-rev(i),1:csz)=m(i)*coefs(drev-rev(i),1:csz)+ &
00311                                            (1-m(i))*coefs(dir,1:csz)
00312                    else ! Fixed node is min
00313                        cfout(drev-rev(i),1:csz)=m(i)*coefs(dir,1:csz)+ &
00314                                            (1-m(i))*coefs(drev-rev(i),1:csz)
00315                    endif
00316                 endif
00317            enddo
00318 
00319        else ! if (nsd==3)
00320 
00321             ! The nodes near the fixed one
00322             do i=1,3
00323                 if (rev(i)>0) then ! Fixed node is max
00324                     cfout(dir+rev(i),1:csz)=m(i)*coefs(dir+rev(i),1:csz)+ &
00325                                         (1-m(i))*coefs(dir,1:csz)
00326                 else ! Fixed node is min
00327                     cfout(dir+rev(i),1:csz)=m(i)*coefs(dir,1:csz)+ &
00328                                         (1-m(i))*coefs(dir+rev(i),1:csz)
00329                 endif         
00330             enddo
00331  
00332             ! The nodes near the new center aka the centers of faces
00333             ! use bilinear interpolation
00334             if (hnds(dh(1))>0) then
00335                 cfout(drev-rev(1),1:csz)=0.0_xyzk; ncnt=ncnt+1
00336                 nc(ncnt)=drev-rev(1); ni(ncnt)=hnds(dh(1))
00337             else
00338                 k=max(0, -rev(1) ) + 1
00339                 cfout(drev-rev(1),1:csz)=m(2)*m(3)*coefs(k,1:csz)+ &
00340                                      m(2)*(1-m(3))*coefs(k+4,1:csz) + &
00341                                      (1-m(2))*m(3)*coefs(k+2,1:csz) + &
00342                                      (1-m(2))*(1-m(3))*coefs(k+6,1:csz)
00343 
00344             endif
00345             
00346             if (hnds(dh(2))>0) then
00347                 cfout(drev-rev(2),1:csz)=0.0_xyzk; ncnt=ncnt+1
00348                 nc(ncnt)=drev-rev(2); ni(ncnt)=hnds(dh(2))
00349             else
00350                 k=max(0, -rev(2) ) + 1
00351                 cfout(drev-rev(2),1:csz)=m(1)*m(3)*coefs(k,1:csz)+ &
00352                                      m(1)*(1-m(3))*coefs(k+4,1:csz) + &
00353                                      (1-m(1))*m(3)*coefs(k+1,1:csz) + &
00354                                      (1-m(1))*(1-m(3))*coefs(k+5,1:csz)
00355            endif
00356 
00357             if (hnds(dh(3))>0) then
00358                 cfout(drev-rev(3),1:csz)=0.0_xyzk; ncnt=ncnt+1
00359                 nc(ncnt)=drev-rev(3); ni(ncnt)=hnds(dh(3))
00360             else
00361                 k=max(0, -rev(3) ) + 1
00362                 cfout(drev-rev(3),1:csz)=m(1)*m(2)*coefs(k,1:csz)+ &
00363                                      m(1)*(1-m(2))*coefs(k+2,1:csz) + &
00364                                      (1-m(1))*m(2)*coefs(k+1,1:csz) + &
00365                                      (1-m(1))*(1-m(2))*coefs(k+3,1:csz)
00366            endif
00367 
00368        endif
00369 
00370 !----------------------------------------------------------------------
00371 ! Following is the previous code with loops and conditionals unrolled for 
00372 ! 2D case. It should be a lot more efficient, in case performance here becomes
00373 ! an issue
00374 
00375        ! NE part
00376 !       if (ct(1)<=pt(1) .and. ct(2)<=pt(2)) then
00377 !          cfout(1,1:csz)=coefs(1,1:csz) 
00378 !          cfout(4,1:csz)=0.0_xyzk; nc(1)=4; ni(1)=ci
00379 
00380 !          if (hnds(2)>0) then
00381 !              cfout(2,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=2; ni(ncnt)=hnds(2)
00382 !          else; cfout(2,1:csz)=m(1)*coefs(2,1:csz)+(1-m(1))*coefs(1,1:csz)
00383 !          endif
00384 
00385 !          if (hnds(1)>0) then
00386 !              cfout(3,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=3; ni(ncnt)=hnds(1)
00387 !          else; cfout(3,1:csz)=m(2)*coefs(3,1:csz)+(1-m(2))*coefs(1,1:csz)
00388 !          endif
00389 
00390        ! SE part
00391 !       else if (ct(1)<=pt(1) .and. ct(2)>pt(2)) then
00392 !          cfout(3,1:csz)=coefs(3,1:csz) 
00393 !          cfout(2,1:csz)=0.0_xyzk; nc(1)=2; ni(1)=ci
00394 
00395 !          if (hnds(4)>0) then
00396 !              cfout(4,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=4; ni(ncnt)=hnds(4)
00397 !          else; cfout(4,1:csz)=m(1)*coefs(4,1:csz)+(1-m(1))*coefs(3,1:csz)
00398 !          endif
00399 
00400 !          if (hnds(1)>0) then
00401 !              cfout(1,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=1; ni(ncnt)=hnds(1)
00402 !          else; cfout(1,1:csz)=m(2)*coefs(3,1:csz)+(1-m(2))*coefs(1,1:csz)
00403 !          endif
00404 
00405        ! NW part
00406 !       else if (ct(1)>pt(1) .and. ct(2)<=pt(2)) then
00407 !          cfout(2,1:csz)=coefs(2,1:csz)
00408 !          cfout(3,1:csz)=0.0_xyzk; nc(1)=3; ni(1)=ci
00409 
00410 !          if (hnds(2)>0) then
00411 !              cfout(1,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=1; ni(ncnt)=hnds(2)
00412 !          else; cfout(1,1:csz)=m(1)*coefs(2,1:csz)+(1-m(1))*coefs(1,1:csz)
00413 !          endif
00414 
00415 !          if (hnds(3)>0) then
00416 !              cfout(4,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=4; ni(ncnt)=hnds(3)
00417 !          else; cfout(4,1:csz)=m(2)*coefs(4,1:csz)+(1-m(2))*coefs(2,1:csz)
00418 !          endif
00419 
00420        ! SW part
00421 !       else if (ct(1)>pt(1) .and. ct(2)>pt(2)) then
00422 !          cfout(4,1:csz)=coefs(4,1:csz)
00423 !          cfout(1,1:csz)=0.0_xyzk; nc(1)=1; ni(1)=ci
00424 
00425 !          if (hnds(4)>0) then
00426 !              cfout(3,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=3; ni(ncnt)=hnds(4)
00427 !          else; cfout(3,1:csz)=m(1)*coefs(4,1:csz)+(1-m(1))*coefs(3,1:csz)
00428 !          endif
00429 
00430 !          if (hnds(3)>0) then
00431 !              cfout(2,1:csz)=0.0_xyzk; ncnt=ncnt+1; nc(ncnt)=2; ni(ncnt)=hnds(3)
00432 !          else; cfout(2,1:csz)=m(2)*coefs(4,1:csz)+(1-m(2))*coefs(2,1:csz)
00433 !          endif
00434 
00435 !       endif
00436 !----------------------------------------------------------------------
00437 
00438        ! Add new nodes to the gaps and thus remove the gaps
00439        do i=ncsz,1,-1
00440          if (all(cfout(:,i)==0.0_xyzk)) then
00441              if (ncnt>0) then ! use one of the new coarse points
00442                  cfout(nc(ncnt),i)=1.0_xyzk; ciout(i)=ni(ncnt); 
00443                  ncnt=ncnt-1
00444              else ! shift one from back to fill the gap
00445                  cfout(:,i)=cfout(:,ncsz)
00446                  ciout(i)=ciout(ncsz)
00447                  ncsz=ncsz-1
00448              endif
00449          endif
00450        enddo
00451 
00452        ! Append all those that didnt fit into gaps
00453        do i=1,ncnt
00454           ncsz=ncsz+1; cfout(:,ncsz)=0.0_xyzk
00455           cfout(nc(i),ncsz)=1.0_xyzk; ciout(ncsz)=ni(i); 
00456        enddo
00457 
00458        if (all(hnds==0)) deallocate(hnds)
00459        
00460 !       write(stream,*) "-------------------------------------"
00461 
00462 !       write(stream,*) "inds:", ACHAR(ciout(1:ncsz)+32)
00463 !       do i=1,2**nsd
00464 !!           if (sum(cfout(i,1:ncsz))<0.9999_xyzk) &
00465 !               write(stream,*) i, ":",cfout(i,1:ncsz)
00466 !       enddo
00467        
00468    end subroutine CalcNextCoefs
00469 
00470    ! Create multipliers for an element bounded by minv/maxv
00471    ! and with corner coefficients of coefs
00472    subroutine BuildMults(minv,maxv,coefs,cinds,csz,nsd,mults)
00473        use RealKind
00474         
00475        real(kind=xyzk),intent(in) :: minv(:),maxv(:)
00476        real(kind=xyzk),intent(in) :: coefs(:,:)
00477        integer, intent(in) :: cinds(:)
00478        integer, intent(in) :: csz
00479        integer, intent(in) :: nsd
00480        real(kind=xyzk), intent(out) :: mults(:,:)
00481 
00482        real(kind=xyzk) :: mbuf(2**nsd), prd
00483        integer :: i, j, k
00484 
00485        mults=0.0_xyzk; prd=1.0_xyzk/product(maxv-minv)
00486 
00487        ! Again, I apologize, but this seemed the most efficient and easy
00488        if (nsd==2) then
00489            ! 1: ne
00490            mbuf(1)=minv(1)*minv(2); mbuf(2)=-minv(2); 
00491            mbuf(3)=-minv(1); mbuf(4)=1.0_xyzk; mbuf=mbuf*prd           
00492            do i=1,csz; mults(:,i)=mults(:,i)+coefs(1,i)*mbuf(:); enddo
00493 
00494            ! 2: nw
00495            mbuf(1)=-maxv(1)*minv(2); mbuf(2)=minv(2)
00496            mbuf(3)=maxv(1); mbuf(4)=-1.0_xyzk; mbuf=mbuf*prd
00497            do i=1,csz; mults(:,i)=mults(:,i)+coefs(2,i)*mbuf(:); enddo
00498 
00499            ! 3: se
00500            mbuf(1)=-minv(1)*maxv(2); mbuf(2)=maxv(2)
00501            mbuf(3)=minv(1); mbuf(4)=-1.0_xyzk; mbuf=mbuf*prd
00502            do i=1,csz; mults(:,i)=mults(:,i)+coefs(3,i)*mbuf(:); enddo
00503 
00504            ! 4: sw
00505            mbuf(1)=maxv(1)*maxv(2); mbuf(2)=-maxv(2); 
00506            mbuf(3)=-maxv(1); mbuf(4)=1.0_xyzk; mbuf=mbuf*prd
00507            do i=1,csz; mults(:,i)=mults(:,i)+coefs(4,i)*mbuf(:); enddo
00508        else ! if (nsd==3)
00509            ! 1: neu
00510            mbuf(5)=-minv(1); mbuf(6)=-minv(2); mbuf(7)=-minv(3)
00511            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=mbuf(1)/mbuf(5:7)
00512            mbuf(8)=1.0_xyzk; mbuf=mbuf*prd
00513            do i=1,csz; mults(:,i)=mults(:,i)+coefs(1,i)*mbuf(:); enddo
00514 
00515            ! 2: nwu
00516            mbuf(5)=maxv(1); mbuf(6)=minv(2); mbuf(7)=minv(3)
00517            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=-mbuf(1)/mbuf(5:7)
00518            mbuf(8)=-1.0_xyzk; mbuf=mbuf*prd
00519            do i=1,csz; mults(:,i)=mults(:,i)+coefs(2,i)*mbuf(:); enddo
00520 
00521            ! 3: seu
00522            mbuf(5)=minv(1); mbuf(6)=maxv(2); mbuf(7)=minv(3)
00523            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=-mbuf(1)/mbuf(5:7)
00524            mbuf(8)=-1.0_xyzk; mbuf=mbuf*prd
00525            do i=1,csz; mults(:,i)=mults(:,i)+coefs(3,i)*mbuf(:); enddo
00526 
00527            ! 4: swu
00528            mbuf(5)=-maxv(1); mbuf(6)=-maxv(2); mbuf(7)=-minv(3)
00529            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=mbuf(1)/mbuf(5:7)
00530            mbuf(8)=1.0_xyzk; mbuf=mbuf*prd
00531            do i=1,csz; mults(:,i)=mults(:,i)+coefs(4,i)*mbuf(:); enddo
00532 
00533            ! 5: ned
00534            mbuf(5)=minv(1); mbuf(6)=minv(2); mbuf(7)=maxv(3)
00535            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=-mbuf(1)/mbuf(5:7)
00536            mbuf(8)=-1.0_xyzk; mbuf=mbuf*prd
00537            do i=1,csz; mults(:,i)=mults(:,i)+coefs(5,i)*mbuf(:); enddo
00538 
00539            ! 6: nwd
00540            mbuf(5)=-maxv(1); mbuf(6)=-minv(2); mbuf(7)=-maxv(3)
00541            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=mbuf(1)/mbuf(5:7)
00542            mbuf(8)=1.0_xyzk; mbuf=mbuf*prd
00543            do i=1,csz; mults(:,i)=mults(:,i)+coefs(6,i)*mbuf(:); enddo
00544 
00545            ! 7: sed
00546            mbuf(5)=-minv(1); mbuf(6)=-maxv(2); mbuf(7)=-maxv(3)
00547            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=mbuf(1)/mbuf(5:7)
00548            mbuf(8)=1.0_xyzk; mbuf=mbuf*prd
00549            do i=1,csz; mults(:,i)=mults(:,i)+coefs(7,i)*mbuf(:); enddo
00550 
00551            ! 8: swd
00552            mbuf(5)=maxv(1); mbuf(6)=maxv(2); mbuf(7)=maxv(3)
00553            mbuf(1)=product(mbuf(5:7)); mbuf(2:4)=-mbuf(1)/mbuf(5:7)
00554            mbuf(8)=-1.0_xyzk; mbuf=mbuf*prd
00555            do i=1,csz; mults(:,i)=mults(:,i)+coefs(8,i)*mbuf(:); enddo
00556 
00557           
00558 !           call DOUG_abort("Trilinear case not implemented yet!")
00559        endif
00560    end subroutine BuildMults
00561 
00562     !! Calculate the multilinear interpolation value for a fine-coarse pair
00563     function MlinInterpolate(pt,nsd,mults) result(res)
00564         use RealKind
00565         real(kind=xyzk), intent(in) :: pt(:)
00566         integer, intent(in) :: nsd
00567         real(kind=xyzk), intent(in) :: mults(:)
00568         real(kind=xyzk) :: res
00569 
00570         ! Bilinear interpolation
00571         if (nsd==2) then
00572            res=mults(1)+ & ! 1
00573                mults(2)*pt(1)+mults(3)*pt(2)+ & ! x/y
00574                mults(4)*pt(1)*pt(2) ! xy
00575         ! Trilinear interpolation
00576         else if (nsd==3) then
00577            res=mults(1)+ & ! 1
00578                mults(2)*pt(1)+mults(3)*pt(2)+mults(4)*pt(3)+ & ! x/y/z
00579                mults(5)*pt(2)*pt(3)+mults(6)*pt(1)*pt(3)+mults(7)*pt(1)*pt(2)+ &
00580                mults(8)*pt(1)*pt(2)*pt(3) ! xyz
00581         endif
00582     end function MLinInterpolate
00583 
00584 
00585    subroutine CalcMlinearInterp(A,C,M)
00586        use RealKind
00587        use Mesh_class
00588        use SpMtx_class
00589        use CoarseGrid_class
00590                 
00591        type(SpMtx),intent(inout) :: A
00592        type(Mesh), intent(in) :: M
00593        type(CoarseGrid), intent(in) :: C
00594 
00595        ! This algorithm calculates practically nothing twice.
00596        ! That can make it hard to follow!
00597        ! The reader should be intimately familiar with the construction of
00598        ! the coarse grid and the layout of its data structures
00599 
00600        ! Initial coarse nodes have their backward edges included
00601        ! (picture assuming positive axes upward and rightward)
00602        ! xxxx
00603        ! i  x
00604        ! i  x
00605        ! iiix       
00606        
00607        ! Refined nodes have the same logic (it is imperative for it to be so!)
00608        ! xxxxx   xxxxx 
00609        ! x x x   x i x
00610        ! xxxxx ; xxiix
00611        ! i x x   x x x
00612        ! iixxx   xxxxx
00613 
00614        ! The element bounds on each level
00615        real(kind=xyzk) :: mins(M%nsd,C%mlvl+1), maxs(M%nsd,C%mlvl+1)
00616 
00617 !       real(kind=xyzk) :: minv(M%nsd), maxv(M%nsd)
00618 
00619        ! NOTICE: 2**M%nsd+C%mlvl+M%nsd is a bound on how many different
00620        ! coarse nodes could be in use for interpolating a single element corner
00621        ! plus (M%nsd+1) spare for simplifying the calculations
00622        
00623        ! coefficcients for coarse nodes for getting each of the element corners
00624        real(kind=xyzk) :: coefs(2**M%nsd,2**M%nsd+C%mlvl+M%nsd,C%mlvl)
00625        ! indices specifing which coarse nodes the coefficients are for
00626        integer :: cinds(2**M%nsd+C%mlvl+M%nsd,C%mlvl)
00627        ! sizes of the previous two arrays
00628        integer :: cszs(C%mlvl)
00629 
00630        ! Have we calculated the final data for interpolation in that dir
00631        logical :: haveData(2**M%nsd)
00632        ! current elements subelement coefficients
00633        real(kind=xyzk) :: curcfs(2**M%nsd,2**M%nsd+C%mlvl+M%nsd,2**M%nsd)
00634        ! its coarse node indices
00635        integer :: curcis(2**M%nsd+C%mlvl+M%nsd,2**M%nsd)
00636        ! and their sizes
00637        integer :: curcszs(2**M%nsd)
00638        ! final values, used directly for interpolation
00639        real(kind=xyzk) :: finals(2**M%nsd,2**M%nsd+C%mlvl+M%nsd,2**M%nsd)
00640        ! Same, but used for where no refinement has occured
00641        real(kind=xyzk) :: mults(2**M%nsd,2**M%nsd)
00642 
00643        real(kind=xyzk),pointer :: ct(:), pt(:)
00644 
00645        integer :: i, j, k, l, o
00646        integer :: d, ci, lvl, cur, pn
00647 
00648        integer :: x, y, z
00649 
00650         cur=1
00651 
00652         do i=1,C%elnum
00653         
00654             mins(:,1)=C%coords(:,C%els(i)%n(1))
00655             maxs(:,1)=C%coords(:,C%els(i)%n(2**M%nsd))
00656 
00657 !            write(stream,*) "/////////////////////////"
00658 !            do j=1,2**M%nsd
00659 !                write (stream,*) j,":",C%coords(:,C%els(i)%n(j))
00660 !            enddo
00661 !            write(stream,*) "/////////////////////////"
00662 
00663             ! Init the level 0 of coefs and cinds
00664             cszs(1)=2**M%nsd
00665             cinds(1:cszs(1),1)=C%els(i)%n(:)
00666             coefs(:,1:cszs(1),1)=0.0_xyzk
00667 
00668             ! As our n-s are chosen in different order
00669             ! One can do a getDir from element center to its corners
00670             !  to check these values, if need be
00671             if (M%nsd==2) then
00672                 coefs(1,4,1)=1.0_xyzk; coefs(2,2,1)=1.0_xyzk
00673                 coefs(3,3,1)=1.0_xyzk; coefs(4,1,1)=1.0_xyzk
00674             else
00675                 coefs(8,1,1)=1.0_xyzk; coefs(4,2,1)=1.0_xyzk
00676                 coefs(6,3,1)=1.0_xyzk; coefs(2,4,1)=1.0_xyzk
00677                 coefs(7,5,1)=1.0_xyzk; coefs(3,6,1)=1.0_xyzk
00678                 coefs(5,7,1)=1.0_xyzk; coefs(1,8,1)=1.0_xyzk
00679             endif
00680 
00681 !       write(stream,*) "/////////////////////////////////////"
00682 
00683 !       write(stream,*) "inds:", ACHAR(cinds(1:cszs(1),1)+32)
00684 !       do k=1,2**M%nsd
00685 !           write(stream,*) k, ":",coefs(k,1:cszs(1),1)
00686 !       enddo
00687  
00688             ! If not refined loop over all the nodes contained in the element
00689             if (C%els(i)%rbeg==-1) then
00690                 ! Generate the final interpolation values
00691                 call BuildMults(mins(:,1),maxs(:,1),coefs(:,:,1),cinds(:,1), &
00692                                 cszs(1),M%nsd,mults)
00693 
00694                 ! Apply them
00695                 do j=C%els(i)%lbeg,C%els(i)%lbeg+C%els(i)%nfs-1
00696                     pn=C%elmap(j) ! The index of the fine node
00697 !                    write(stream,*) pn, "> IndsInitial:",cinds(:,0)," cur:",cur
00698 
00699                     ! Calculate the values
00700                     l=1
00701                     do k=cur,cur+cszs(1)-1
00702 
00703                         A%val( k )= &
00704                              MlinInterpolate(M%lcoords(:,pn),M%nsd,mults(:,l))
00705                         A%indi(k)=pn; A%indj(k)=cinds(l,1)
00706                         l=l+1
00707                     enddo
00708                 
00709                     ! Advance the pointer
00710                     cur=cur+cszs(1)
00711                 enddo
00712 
00713             else 
00714 
00715             ! Otherwise loop through the refinements
00716             j=C%els(i)%rbeg
00717             do while (j/=-1)
00718 
00719                 !write(stream,*),"going to ",C%refels(j)%level
00720 
00721                 ! Calculate the previous level data correctly
00722                 if (C%refels(j)%parent>0) then ! not first
00723                     ci=C%refels(C%refels(j)%parent)%node
00724                     ct=>C%coords(:,ci); lvl=C%refels(j)%level-1
00725                     pt=>C%coords(:,C%refels(j)%node)
00726 !                    write(stream,*) "LEVEL:",C%refels(C%refels(j)%parent)%level
00727                     call CalcNextCoefs(ct,ci,&
00728                              C%refels(C%refels(j)%parent)%hnds,&
00729                              mins(:,lvl),maxs(:,lvl),&
00730                              coefs(:,:,lvl),cinds(:,lvl),cszs(lvl),M%nsd,pt, &
00731                              coefs(:,:,lvl+1),cinds(:,lvl+1),cszs(lvl+1))
00732 
00733                     ! Adjust the bounds
00734                     do l=1,M%nsd
00735                         if (ct(l)<=pt(l)) then
00736                             mins(l,lvl+1)=ct(l); maxs(l,lvl+1)=maxs(l,lvl)
00737                         else
00738                             mins(l,lvl+1)=mins(l,lvl); maxs(l,lvl+1)=ct(l)
00739                         endif
00740                     enddo
00741 !                    write(stream,*) "Mins: ",mins(:,lvl+1)
00742 !                    write(stream,*) "Maxs: ",maxs(:,lvl+1)
00743 
00744 ! For testing getRefBounds
00745 !                    call getRefBounds(j,C,M%nsd,minv,maxv)
00746 !                    if (any(minv/=mins(:,lvl+1)) .or. any(maxv/=maxs(:,lvl+1))) then
00747 !                        write (stream,*) "XX*****************************"
00748 !                        write (stream,*) minv
00749 !                        write (stream,*) mins(:,lvl+1)
00750 !                        write (stream,*) maxv
00751 !                        write (stream,*) maxs(:,lvl+1)
00752 !                    else 
00753 !                        write (stream,*) "Bounds OK"
00754 !                    endif
00755                         
00756                 endif   
00757 
00758                 ci=C%refels(j)%node; ct=>C%coords(:,ci)
00759                 lvl=C%refels(j)%level
00760 
00761                 ! Check if we actually need to do anything here
00762                 if (C%refels(j)%lbeg<=C%refels(j)%lend) then
00763                     haveData=.false.
00764        
00765                     ! Go through all the fine nodes that are in this element
00766                     do k=C%refels(j)%lbeg,C%refels(j)%lend
00767                         pn=C%elmap(k); pt=>M%lcoords(:,pn)
00768 
00769                         ! Find the direction
00770                         d=getDir(pt-ct,M%nsd)
00771 
00772                         ! Calculate the values if nessecary
00773                        if (.not.haveData(d)) then
00774                            call CalcNextCoefs(ct,ci,&
00775                                 C%refels(j)%hnds,&
00776                                 mins(:,lvl),maxs(:,lvl),&
00777                                 coefs(:,:,lvl),cinds(:,lvl),cszs(lvl), &
00778                                 M%nsd,pt,curcfs(:,:,d),curcis(:,d), curcszs(d))
00779 
00780                            ! Adjust the bounds
00781                            do l=1,M%nsd
00782                                if (ct(l)<=pt(l)) then
00783                                   mins(l,lvl+1)=ct(l); maxs(l,lvl+1)=maxs(l,lvl)
00784                                else
00785                                   mins(l,lvl+1)=mins(l,lvl); maxs(l,lvl+1)=ct(l)
00786                                endif
00787                            enddo
00788 
00789 !                            write(stream,*) "Mins: ",mins(:,lvl+1)
00790 !                            write(stream,*) "Maxs: ",maxs(:,lvl+1)
00791 
00792                            call BuildMults( mins(:,lvl+1),maxs(:,lvl+1), &
00793                                curcfs(:,:,d),curcis(:,d),curcszs(d), &
00794                                M%nsd,finals(:,:,d))
00795 
00796                         
00797                            haveData(d)=.true.
00798                        endif
00799 
00800                        ! Use the values for interpolation 
00801                        l=1
00802                        do o=cur,cur+curcszs(d)-1
00803                           A%val( o )= &
00804                             MlinInterpolate(M%lcoords(:,pn),M%nsd,finals(:,l,d))
00805 
00806                           A%indi(o)=pn; A%indj(o)=curcis(l,d)
00807                           l=l+1
00808                        enddo
00809                 
00810                        ! Advance the pointer
00811                        cur=cur+curcszs(d)
00812                     enddo
00813                 endif
00814 
00815                 ! And go to the next refined element
00816                 j=C%refels(j)%next
00817 
00818             enddo
00819             endif
00820         enddo
00821 
00822         ! Mark down how much we actually used 
00823         A%nnz=cur-1
00824 
00825    end subroutine
00826 
00827 end module GeomInterp