|
DOUG 0.2
|
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 00043 module CoarseCreateRestrict 00044 use RealKind 00045 00046 #include<doug_config.h> 00047 00048 #ifdef D_COMPLEX 00049 #define float complex 00050 #else 00051 #define float real 00052 #endif 00053 ! For now 00054 real(kind=xyzk), parameter :: invdistpow=0.5_xyzk 00055 00056 contains 00057 00058 !! Create the Restriction matrix 00059 subroutine CreateRestrict(C,M,R) 00060 use CoarseGrid_class 00061 use GeomInterp 00062 use SpMtx_class 00063 use CoarseMtx_mod 00064 use globals 00065 00066 implicit none 00067 00068 !! Coarse Grid whose structure to use 00069 type(CoarseGrid), intent(inout) :: C 00070 !! Fine Mesh for which to build 00071 type(Mesh), intent(in) :: M 00072 !! Restriction matrix to create 00073 type(SpMtx), intent(out) :: R 00074 00075 00076 ! We can add a choice here to do different things 00077 00078 select case(sctls%interpolation_type) 00079 00080 case (COARSE_INTERPOLATION_INVDIST) 00081 call CreatePathRestrict(C,M,R,genNoData,& 00082 getInvDistVals,getSizeOne) 00083 00084 case (COARSE_INTERPOLATION_KRIGING) 00085 call CreatePathRestrict(C,M,R,genKrigingData,& 00086 getKrigingVals,getKrigingSize) 00087 00088 case (COARSE_INTERPOLATION_MULTLIN) 00089 !call CreatePathRestrict(C,M,Restrict,genMultiLinearData,& 00090 ! getMultiLinearVals,getMultiLinearSize) 00091 call CreateGeneralRestrict(C,M,R,CalcMlinearInterp) 00092 00093 end select 00094 00095 end subroutine 00096 00097 !! Calculates the node prolongation matrix taking the coarse element path 00098 !! down to the fine node. Later rebuilds it into a freedom restrict mat. 00099 subroutine CreatePathRestrict(C,M,R,gendata,getvals,getsize) 00100 use RealKind 00101 use CoarseGrid_class 00102 use SpMtx_class 00103 use SpMtx_util 00104 use globals 00105 00106 implicit none 00107 00108 !! Coarse Grid whose structure to use 00109 type(CoarseGrid), intent(inout) :: C 00110 !! Fine Mesh for which to build 00111 type(Mesh), intent(in) :: M 00112 !! The restriction matrix to be created 00113 type(SpMtx), intent(out) :: R 00114 00115 !! gendata is used to create aux data of size getsize 00116 !! for interpolation from points in pts 00117 interface 00118 subroutine gendata(cpts,csz,nsd,routp,ioutp) 00119 use RealKind 00120 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz) 00121 integer :: csz, nsd 00122 real(kind=xyzk), intent(out) :: routp(:) ! routp(undet.) 00123 integer, intent(out) :: ioutp(:) ! ioutp(undet.) 00124 end subroutine gendata 00125 end interface 00126 00127 !! getvals is used to create the prolongation data for a fine point 00128 interface 00129 subroutine getvals(cpts,csz,nsd,rinp,iinp,pt,outp) 00130 use RealKind 00131 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz) 00132 integer :: csz, nsd 00133 real(kind=xyzk), intent(in) :: rinp(:) ! rinp(undet.) 00134 integer, intent(in) :: iinp(:) ! iinp(undet.) 00135 real(kind=xyzk), intent(in) :: pt(:) ! pt(nsd) 00136 float(kind=rk), intent(out) :: outp(:) ! outp(csz) 00137 end subroutine getvals 00138 end interface 00139 00140 !! get the max size of aux data needed for ptcnt points 00141 interface 00142 subroutine getsize(ptcnt,nsd,rsz,isz) 00143 integer, intent(in) :: ptcnt 00144 integer, intent(in) :: nsd 00145 integer, intent(out) :: rsz, isz 00146 end subroutine getsize 00147 end interface 00148 00149 00150 type(SpMtx) :: NP ! node prolongation matrix - restrict is transposed 00151 00152 integer :: cnt 00153 integer :: i, j, k, pn, nd, f 00154 integer :: tnsd 00155 real(kind=xyzk) :: pts(M%nsd,2**M%nsd+C%mlvl) 00156 integer :: pinds(2**M%nsd+C%mlvl) 00157 real(kind=xyzk),pointer :: pt(:) 00158 real(kind=xyzk),pointer :: raux(:) 00159 integer,pointer :: iaux(:) 00160 00161 ! For the reverse of freemap 00162 integer :: flist(C%nlfc) 00163 integer :: fbeg(C%nct+1) 00164 integer :: faux(C%nct) 00165 00166 tnsd=2**M%nsd 00167 00168 !**************************************************** 00169 ! Create the Node Prolongation matrix structure 00170 !**************************************************** 00171 00172 ! Calculate nnz 00173 00174 cnt=0 00175 00176 ! Increment by initial grid 00177 do i=1, C%elnum 00178 cnt=cnt+tnsd*C%els(i)%nfs 00179 enddo 00180 00181 ! Increment by refinements 00182 do i=1, C%refnum 00183 cnt = cnt + & 00184 C%refels(i)%level*(C%refels(i)%lend-C%refels(i)%lbeg+1) 00185 enddo 00186 00187 ! Allocate the node prolong matrix 00188 NP=SpMtx_newInit(cnt,nrows=M%lnnode,ncols=C%nct) 00189 allocate(NP%M_bound(M%lnnode+1)) 00190 00191 ! Create NP%M_bound - done in 2 phases 00192 NP%M_bound=0 00193 00194 ! Vector in coarse non-refined elements 00195 do i=1, C%elnum 00196 if (C%els(i)%rbeg==-1) then 00197 do j=C%els(i)%lbeg,C%els(i)%lbeg+C%els(i)%nfs-1 00198 NP%M_bound(C%elmap(j)+1)=tnsd 00199 enddo 00200 endif 00201 enddo 00202 00203 ! Vector in refined elements 00204 do i=1, C%refnum 00205 cnt=tnsd+C%refels(i)%level 00206 do j=C%refels(i)%lbeg,C%refels(i)%lend 00207 NP%M_bound(C%elmap(j)+1)=cnt 00208 enddo 00209 enddo 00210 00211 ! And add the counts together to get the final M_bound 00212 NP%M_bound(1)=1 00213 do i=1,M%lnnode 00214 NP%M_bound(i+1)=NP%M_bound(i)+NP%M_bound(i+1) 00215 enddo 00216 00217 !**************************************************** 00218 ! Create the prolongation data for NP 00219 !**************************************************** 00220 00221 ! Init the aux buffer to max possible size 00222 call getsize(C%mlvl+tnsd,M%nsd,i, k) 00223 allocate(raux(i),iaux(k)) 00224 00225 do i=1,C%elnum 00226 00227 ! Gather the elements corner points to an array 00228 do k=1,2**M%nsd 00229 pts(:,k)=C%coords(:,C%els(i)%n(k)) 00230 enddo 00231 00232 ! If not refined loop over all the nodes contained in the element 00233 if (C%els(i)%rbeg==-1) then 00234 ! Generate the auxilliary data 00235 call gendata(pts,tnsd,M%nsd,raux,iaux) 00236 pinds(1:tnsd)=C%els(i)%n 00237 00238 do j=C%els(i)%lbeg,C%els(i)%lbeg+C%els(i)%nfs-1 00239 pn=C%elmap(j) ! The index of the fine node 00240 00241 ! Evaluate the functions in the point 00242 call getvals(pts,tnsd,M%nsd,raux,iaux,M%lcoords(:,pn), & 00243 NP%val(NP%M_bound(pn):NP%M_bound(pn+1)-1)) 00244 00245 ! write(stream,*) "NP",pn,M%lnnode,NP%M_bound(pn),NP%nnz 00246 00247 ! Set appropriate indi and indj 00248 NP%indi(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pn 00249 NP%indj(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pinds 00250 enddo 00251 else ! Otherwise loop through the refinements 00252 pinds(1:tnsd)=C%els(i)%n 00253 j=C%els(i)%rbeg 00254 do while (j/=-1) 00255 ! Set the node of the element into pts 00256 pts(:,tnsd+C%refels(j)%level)=C%coords(:,C%refels(j)%node) 00257 pinds(tnsd+C%refels(j)%level)=C%refels(j)%node 00258 00259 ! Check if we actually need to do anything here 00260 if (C%refels(j)%lbeg<=C%refels(j)%lend) then 00261 call gendata(pts,tnsd+C%refels(j)%level,M%nsd,raux,iaux) 00262 00263 ! Go through all the fine nodes that are in this element 00264 do k=C%refels(j)%lbeg,C%refels(j)%lend 00265 pn=C%elmap(k) 00266 00267 ! Evaluate the functions in the point 00268 call getvals(pts,tnsd+C%refels(j)%level,& 00269 M%nsd,raux,iaux,M%lcoords(:,pn), & 00270 NP%val(NP%M_bound(pn):NP%M_bound(pn+1)-1)) 00271 00272 ! Set appropriate indi and indj 00273 NP%indi(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pn 00274 NP%indj(NP%M_bound(pn):NP%M_bound(pn+1)-1)=pinds 00275 enddo 00276 endif 00277 00278 ! And go to the next refined element 00279 j=C%refels(j)%next 00280 enddo 00281 endif 00282 enddo 00283 00284 deallocate(raux,iaux) 00285 00286 00287 !**************************************************** 00288 ! Create the reverse of freemap 00289 !**************************************************** 00290 00291 ! Find the counts of freedoms associated with nodes 00292 faux=0; 00293 do i=1,C%nlfc 00294 faux(C%cfreemap(i))=faux(C%cfreemap(i))+1 00295 enddo 00296 00297 ! Create fbegs 00298 fbeg(1)=1; 00299 do i=1,C%nct 00300 fbeg(i+1)=faux(i)+fbeg(i) 00301 enddo 00302 00303 ! Create flist 00304 faux=fbeg(1:C%nlfc) 00305 do i=1,C%nlfc 00306 nd=C%cfreemap(i) 00307 flist(faux(nd))=i 00308 faux(nd)=faux(nd)+1 00309 enddo 00310 00311 !**************************************************** 00312 ! Create R based on NP 00313 !**************************************************** 00314 00315 ! Calculate the upper bound for nnz 00316 cnt=0 00317 do i=1,NP%nnz 00318 cnt=cnt+fbeg(NP%indj(i)+1)-fbeg(NP%indj(i)) 00319 enddo 00320 00321 ! Allocate the matrix 00322 R=SpMtx_newInit(cnt,nrows=C%nlfc,ncols=M%nlf) 00323 allocate(R%M_bound(M%nlf+1)) 00324 00325 ! Fill it 00326 f=1; 00327 do i=1,M%nlf 00328 nd=M%lfreemap(i) 00329 R%M_bound(i)=f 00330 do j=NP%M_bound(nd),NP%M_bound(nd+1)-1 00331 do k=fbeg(NP%indj(j)),fbeg(NP%indj(j)+1)-1 00332 if (NP%val(j)>eps) then 00333 R%val(f)=NP%val(j) 00334 R%indi(f)=flist(k) 00335 R%indj(f)=i 00336 f=f+1 00337 endif 00338 enddo 00339 enddo 00340 enddo 00341 R%M_bound(M%nlf+1)=f 00342 00343 ! Set the matrix to be row-sorted format 00344 R%arrange_type=D_SpMtx_ARRNG_COLS 00345 00346 ! Destroy the node version 00347 call SpMtx_Destroy(NP) 00348 00349 ! Give back the unneeded memory in P 00350 if (R%nnz>f-1) call SpMtx_resize(R,f-1) 00351 00352 if (sctls%verbose>3) & 00353 write (stream,*) "Restriction matrix has ",f-1," elements" 00354 00355 end subroutine CreatePathRestrict 00356 00357 !************************************************ 00358 ! Inverse Distances Interpolation 00359 !************************************************ 00360 00361 !! A function to pass to CreateProlong to return 1 00362 subroutine getSizeOne(ptcnt, nsd, rsz, isz) 00363 integer, intent(in) :: ptcnt, nsd 00364 integer, intent(out) :: rsz, isz 00365 rsz=1; isz=1; 00366 end subroutine getSizeOne 00367 00368 !! A function to pass to CreateProlong to create no real data 00369 subroutine genNoData(cpts,csz,nsd,routp,ioutp) 00370 use RealKind 00371 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz) 00372 integer :: csz, nsd 00373 real(kind=xyzk), intent(out) :: routp(:) ! outp(undet.) 00374 integer, intent(out) :: ioutp(:) 00375 routp(1)=0.0_rk; ioutp(1)=0 00376 end subroutine genNoData 00377 00378 !! Get the interpolation values by using inverse distances method 00379 subroutine getInvDistVals(cpts,csz,nsd,rinp,iinp,pt,outp) 00380 use RealKind 00381 use globals, only: eps 00382 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz) 00383 integer :: csz, nsd 00384 real(kind=xyzk), intent(in) :: rinp(:) ! rinp(1) 00385 integer, intent(in) :: iinp(:) ! iinp(1) 00386 real(kind=xyzk), intent(in) :: pt(:) ! pt(nsd) 00387 float(kind=rk), intent(out) :: outp(:) ! outp(csz) 00388 00389 integer :: i 00390 00391 ! Get the distances 00392 do i=1,csz 00393 outp(i)=(dot_product((cpts(:,i)-pt),(cpts(:,i)-pt)))**invdistpow 00394 00395 if (outp(i)<=eps) then 00396 outp(1:csz)=0.0_rk; 00397 outp(i)=1.0_rk; 00398 exit 00399 else 00400 outp(i)=1.0_rk/outp(i) 00401 endif 00402 enddo 00403 00404 ! Normalize to preserve constant 00405 outp(1:csz)=outp(1:csz)/sum(outp(1:csz)) 00406 00407 end subroutine getInvDistVals 00408 00409 !************************************************ 00410 ! Kriging Interpolation 00411 !************************************************ 00412 00413 !! A function to pass to CreateProlong for Kriging 00414 subroutine getKrigingSize(ptcnt, nsd, rsz, isz) 00415 integer, intent(in) :: ptcnt, nsd 00416 integer, intent(out) :: rsz, isz 00417 rsz=(ptcnt+1)**2 ! Matrix 00418 isz=ptcnt*2 ! Pivot data 00419 end subroutine getKrigingSize 00420 00421 !! A function to pass to CreateProlong to create the Kriging aux data 00422 subroutine genKrigingData(cpts,csz,nsd,routp,ioutp) 00423 use RealKind 00424 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz) 00425 integer :: csz, nsd 00426 real(kind=xyzk), intent(out) :: routp(:) ! routp((csz+1)**2) 00427 integer, intent(out) :: ioutp(:) ! ioutp(2*csz) 00428 00429 real(kind=xyzk) :: mat(csz+1,csz+1) 00430 integer :: i, k 00431 00432 ! Create the matrix 00433 do i=1,csz 00434 mat(csz,i)=1.0_xyzk ! Last row 00435 mat(i,csz)=1.0_xyzk ! Last column 00436 mat(i,i)=0.0_xyzk ! Diagonal 00437 00438 ! Other elements 00439 do k=i+1,csz 00440 ! Set the distance 00441 mat(i,k)= & 00442 (dot_product((cpts(:,i)-cpts(:,k)),(cpts(:,i)-cpts(:,k))))**invdistpow 00443 ! Copy it 00444 mat(k,i)=mat(i,k) 00445 enddo 00446 enddo 00447 00448 ! LUP-factorize mat 00449 00450 !!!! call DGETRF(csz+1,csz+1,mat,csz+1,ioutp,i) 00451 00452 ! if (i/=0) ERROR 00453 00454 ! Pack mat into outp 00455 routp(1:(csz+1)**2)=reshape(mat,(/ (csz+1)**2 /) ) 00456 00457 end subroutine genKrigingData 00458 00459 !! Get the interpolation values by using Kriging method 00460 subroutine getKrigingVals(cpts,csz,nsd,rinp,iinp,pt,outp) 00461 use RealKind 00462 real(kind=xyzk), intent(in) :: cpts(:,:) ! pts(nsd,csz) 00463 integer :: csz, nsd 00464 float(kind=rk), intent(in) :: rinp(:) ! rinp((csz+1)**2) 00465 integer, intent(in) :: iinp(:) ! iinp(2*csz) 00466 real(kind=xyzk), intent(in) :: pt(:) ! pt(nsd) 00467 float(kind=rk), intent(out) :: outp(:) ! outp(csz) 00468 00469 float(kind=rk) :: vec(csz+1,1) 00470 integer :: i 00471 00472 ! Calculate the distances vector for this point 00473 do i=1,csz 00474 vec(i,1)=(dot_product((cpts(:,i)-pt),(cpts(:,i)-pt)))**invdistpow 00475 enddo 00476 vec(csz+1,1)=1.0_xyzk 00477 00478 ! Solve it 00479 !!!!call DGETRS('N', csz+1, 1, reshape(rinp,(/ csz+1, csz+1 /)),csz+1, & 00480 !!!! iinp, vec, csz+1, i) 00481 00482 ! if (i/=0) ERROR 00483 00484 ! Copy the data over 00485 outp(1:csz)=vec(1:csz,1) 00486 00487 end subroutine getKrigingVals 00488 00489 00490 !! Strip the restriction matrix 00491 !! useful for multiple processor case, where 00492 !! the stripped variant is used for vector restrict/interpolate 00493 subroutine stripRestrict(M, R) 00494 use Mesh_class 00495 use SpMtx_class 00496 use globals, only : stream 00497 00498 type(Mesh), intent(in) :: M 00499 type(SpMtx), intent(inout) :: R 00500 ! type(SpMtx), intent(out) :: Ra 00501 00502 integer :: i,j,cnt 00503 integer :: unq(M%nlf) 00504 00505 ! Calculate unq (has 1 if that node is "unique" 00506 ! a.k.a. this is the lowest ranked process which has it) 00507 unq=0; cnt=0 00508 do i=1,M%nell 00509 if (M%eptnmap(i)<=myrank) then ! lower ranked 00510 do j=1,M%nfrelt(i) 00511 if (M%gl_fmap(M%mhead(j,i))/=0) & ! non-unique 00512 unq(M%gl_fmap(M%mhead(j,i)))=-1 00513 enddo 00514 elseif (M%eptnmap(i)==myrank+1) then ! my nodes 00515 do j=1,M%nfrelt(i) 00516 if (unq(M%gl_fmap(M%mhead(j,i)))==0) & ! untagged 00517 unq(M%gl_fmap(M%mhead(j,i)))=1 00518 enddo 00519 endif 00520 enddo 00521 00522 ! Count the nnz in the stripped restrict matrix 00523 ! cnt=0 00524 ! if (R%arrange_type==D_SpMtx_ARRNG_COLS) then ! can use M_bound 00525 ! do i=1,R%ncols 00526 ! if (unq(i)==1) cnt=cnt+R%M_bound(i+1)-R%M_bound(i) 00527 ! enddo 00528 ! else ! count the matrix elements one by one 00529 ! do i=1,R%nnz 00530 ! if (unq(R%indj(i))==1) cnt=cnt+1 00531 ! enddo 00532 ! endif 00533 00534 ! ! Initalize Ra 00535 ! Ra=SpMtx_newInit(nnz=cnt,nrows=R%nrows,ncols=R%ncols) 00536 00537 ! if we have something to actually do 00538 ! if (cnt/=0) then 00539 ! Fille Ra with its elements 00540 j=1 00541 do i=1,R%nnz 00542 if (unq(R%indj(i))==1) then ! it stays 00543 R%indi(j)=R%indi(i) 00544 R%indj(j)=R%indj(i) 00545 R%val(j)=R%val(i) 00546 j=j+1 00547 endif 00548 enddo 00549 ! endif 00550 00551 call SpMtx_resize(R,j-1) 00552 R%arrange_type=D_SpMtx_ARRNG_NO 00553 deallocate(R%M_bound) 00554 00555 if (sctls%verbose>3) & 00556 write (stream,*) "Stripped restriction matrix has ",j-1," elements" 00557 00558 end subroutine stripRestrict 00559 00560 end module CoarseCreateRestrict
1.7.3-20110217