DOUG 0.2

ElemMtxs_assemble.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 !-------------------------------------------------------
00025 module ElemMtxs_assemble
00026 
00027   use SpMtx_class
00028   use SpMtx_permutation
00029   use ElemMtxs_base
00030   use IdxMap_class
00031   use Mesh_class
00032   use RealKind
00033   use globals
00034 
00035   implicit none
00036 
00037 #include<doug_config.h>
00038 
00039 #ifdef D_COMPLEX
00040 #define float complex
00041 #else
00042 #define float real
00043 #endif
00044 
00045   !---------------------------------------------
00049   type ElemMtxsAssembleBlock
00051      integer :: val_count
00052 
00053      float(kind=rk), dimension(:), pointer :: val
00055      type(IdxMap), dimension(:), pointer :: idx_map
00056   end type ElemMtxsAssembleBlock
00057 
00058 
00059   !---------------------------------------------
00063   type ElemMtxsAssembleContext
00065      type(ElemMtxsAssembleBlock), dimension(2, 2) :: blocks
00066 
00067      float(kind=rk), dimension(:), pointer     :: rhs
00069      integer, dimension(:), pointer :: inv_perm_map
00070   end type ElemMtxsAssembleContext
00071 
00072   !private :: &
00073   !     ElemMtxsAssembleBlock
00074 
00075 
00076 contains
00077 
00078 
00079   !----------------------------------------------------------------
00082   function ElemMtxsAssembleContext_newInit(Msh) result(AC)
00083     implicit none
00084 
00085     type(Mesh),        intent(in) :: Msh
00086     type(ElemMtxsAssembleContext),target :: AC
00087 
00088     integer :: i, m, n
00089     integer :: iin, iintf, nfintf
00090     type(ElemMtxsAssembleBlock), pointer :: AB
00091 
00092     ! Calculate inverse permuation map
00093     allocate(AC%inv_perm_map(Msh%nlf))
00094     nfintf = sum(Msh%inner_interf_fmask) ! number of interface freedoms
00095     iin = nfintf
00096     iintf = 0
00097     do i = 1,Msh%nlf
00098        if (Msh%inner_interf_fmask(i) == D_FREEDOM_INNER) then
00099           iin = iin + 1
00100           AC%inv_perm_map(i) = iin
00101        else
00102           iintf = iintf + 1
00103           AC%inv_perm_map(i) = iintf
00104        end if
00105     end do
00106     
00107     ! Allocate space for RHS
00108     allocate(AC%rhs(Msh%nlf))
00109     AC%rhs = 0.0_rk
00110 
00111     ! Initialize blocks
00112     do n = 1, 2
00113        do m = 1, 2
00114           AB => AC%blocks(n,m)
00115           AB%val => NULL()
00116           AB%val_count = 0
00117           allocate(AB%idx_map(Msh%nlf)) !! TODO: make this global for context?
00118           do i = 1,Msh%nlf
00119              AB%idx_map(i) = IdxMap_New() !! TODO: make this global for context?
00120           end do
00121        end do
00122     end do
00123   end function ElemMtxsAssembleContext_newInit
00124 
00125 
00126   !----------------------------------------------------------------
00129   subroutine ElemMtxsAssembleContext_Destroy(AC)
00130     implicit none
00131 
00132     type(ElemMtxsAssembleContext), intent(in out),target :: AC
00133 
00134     integer :: i, m, n
00135     type(ElemMtxsAssembleBlock), pointer :: AB
00136 
00137     ! Free all blocks
00138     do n = 1, 2
00139        do m = 1, 2
00140           AB => AC%blocks(n,m)
00141           do i = 1,size(AB%idx_map)
00142              call IdxMap_Destroy(AB%idx_map(i))
00143           end do
00144           deallocate(AB%idx_map)
00145           if (associated(AB%val)) deallocate(AB%val)
00146        end do
00147     end do
00148 
00149     ! Free other resources
00150     deallocate(AC%rhs)
00151     deallocate(AC%inv_perm_map)
00152   end subroutine ElemMtxsAssembleContext_Destroy
00153 
00154 
00155   !----------------------------------------------------------------
00158   subroutine ElemMtxsAssembleContext_addChunk(AC, E, Msh)
00159     implicit none
00160 
00161     type(ElemMtxsAssembleContext), intent(in out),target :: AC  ! Assemble context
00162     type(ElemMtxsChunk),        intent(in)     :: E   ! Element matrices and RHSs
00163     type(Mesh),                 intent(in)     :: Msh ! Mesh
00164 
00165     integer :: m, n, i, j, indi, indj, idx, val_size
00166     integer :: ge, le, gf, lf, colgf, collf
00167     float(kind=rk), dimension(:), pointer :: val_temp
00168     type(ElemMtxsAssembleBlock), pointer :: AB
00169 
00170     do le = 1,E%nell
00171        ge = E%lg_emap(le)
00172        do i = 1,Msh%nfrelt(ge)
00173           gf = Msh%mhead(i,ge)
00174           lf = Msh%gl_fmap(gf)
00175           if (lf > 0) then ! if the node is mine at all?
00176              ! Assemble elem part
00177              n = 1
00178              if (Msh%inner_interf_fmask(lf) /= D_FREEDOM_INTERF) n = 2
00179              do j = 1,Msh%nfrelt(ge)
00180                 if (E%elem(i, j, le) == 0.0_rk) cycle
00181                 colgf = Msh%mhead(j,ge)
00182                 collf = Msh%gl_fmap(colgf)
00183                 if (collf > 0) then ! if the node is mine at all?
00184                    m = 1
00185                    if (Msh%inner_interf_fmask(collf) /= D_FREEDOM_INTERF) m = 2
00186                    AB => AC%blocks(n,m)
00187                    indi = AC%inv_perm_map(lf)
00188                    indj = AC%inv_perm_map(collf)
00189                    idx = IdxMap_Lookup(AB%idx_map(indi), indj)
00190                    if (idx == -1) then
00191                       AB%val_count = AB%val_count + 1
00192                       idx = AB%val_count
00193                       val_size = 0
00194                       if (associated(AB%val)) val_size = size(AB%val)
00195                       if (val_size < idx) then
00196                          allocate(val_temp(idx * 4 / 3 + 16))
00197                          val_temp = 0.0_rk
00198                          if (associated(AB%val)) then
00199                             val_temp(1:val_size) = AB%val(1:val_size)
00200                             deallocate(AB%val)
00201                          end if
00202                          AB%val => val_temp
00203                       end if
00204                       call IdxMap_Insert(AB%idx_map(indi), indj, idx)
00205                    end if
00206                    AB%val(idx) = AB%val(idx) + E%elem(i, j, le)
00207                 end if
00208              end do
00209              
00210              ! Assemble RHS part
00211              indi = AC%inv_perm_map(lf)
00212              AC%rhs(indi) = AC%rhs(indi) + E%elemrhs(i, le)
00213           end if
00214        end do
00215     end do
00216   end subroutine ElemMtxsAssembleContext_addChunk
00217 
00218 
00219   !----------------------------------------------------------------
00229   subroutine ElemMtxsAssembleContext_extractSpMtx(AC, A, Msh)
00230     use globals
00231     implicit none
00232   
00233     type(ElemMtxsAssembleContext), intent(in),target  :: AC
00234     type(SpMtx),                intent(out) :: A !< system matrix, should be uninitialized before calling
00235     type(Mesh),                 intent(in)  :: Msh
00236 
00237     integer :: i, j, k, n, m
00238     integer :: nnz, indi, indj, idx, bbe
00239     type(ElemMtxsAssembleBlock), pointer :: AB
00240 
00241     ! Calculate number of non-zero elements
00242     nnz = 0
00243     do m=1,2
00244        do n=1,2
00245           nnz = nnz + AC%blocks(n,m)%val_count
00246        end do
00247     end do
00248     
00249     ! Create sparse matrix object
00250     A = SpMtx_newInit(nnz, sctls%number_of_blocks, Msh%nlf, &
00251          symmstruct=sctls%symmstruct, symmnumeric=sctls%symmnumeric)
00252     call SpMtx_buildPermMap(A, Msh)
00253     
00254     ! Fill sparse matrix object
00255     bbe = 1
00256     do m = 1,2
00257        do n = 1,2
00258           AB => AC%blocks(n,m)
00259           A%mtx_bbs(n,m) = bbe
00260           do i = 1,Msh%nlf
00261              do k = 1,IdxMap_Size(AB%idx_map(i))
00262                 j = IdxMap_Key(AB%idx_map(i), k)
00263                 idx = IdxMap_Lookup(AB%idx_map(i), j)
00264                 A%val(bbe)  = AB%val(idx)
00265                 A%indi(bbe) = j ! this flipping is needed to get identical behaviour to old DOUG
00266                 A%indj(bbe) = i
00267                 bbe = bbe + 1
00268              end do
00269           end do
00270           A%mtx_bbe(n,m) = bbe-1
00271        end do
00272     end do
00273 
00274     call SpMtx_setMtxInnerBound(A, A%mtx_bbs(2,2))
00275   end subroutine ElemMtxsAssembleContext_extractSpMtx
00276 
00277 
00278   !----------------------------------------------------------------
00281   subroutine ElemMtxsAssembleContext_extractVect(AC, b, Msh)
00282     implicit none
00283 
00284     type(ElemMtxsAssembleContext),   intent(in)  :: AC
00285     float(kind=rk), dimension(:), intent(out) :: b
00286     type(Mesh),                   intent(in)  :: Msh
00287 
00288     b = AC%rhs
00289   end subroutine ElemMtxsAssembleContext_extractVect
00290 
00291 end module ElemMtxs_assemble