From 081df995153f249fc86c04972692d43b2f7450dd Mon Sep 17 00:00:00 2001 From: miguel Date: Thu, 29 Mar 2001 16:02:03 +0000 Subject: Removed the original file that handled extrapolation directly. There are now three files for three purposes: findboundary takes a mask with only 0s and 1s and marks the boundary with 0.5s. findnormals takes a mask where the boundary has been marked and returns the normal directions to the mask for these locations. extrapolate takes a mask the normals and extrapolates a grid function. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/LegoExcision/trunk@6 f75ba9e5-694f-0410-ac2c-87ea7ce7132b --- src/extrapolate.F90 | 62 +++++++++++ src/findboundary.F90 | 57 ++++++++++ src/findnormals.F90 | 290 +++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 2 +- src/maskvalues.h | 8 ++ 5 files changed, 418 insertions(+), 1 deletion(-) create mode 100644 src/extrapolate.F90 create mode 100644 src/findboundary.F90 create mode 100644 src/findnormals.F90 create mode 100644 src/maskvalues.h diff --git a/src/extrapolate.F90 b/src/extrapolate.F90 new file mode 100644 index 0000000..b9b0b99 --- /dev/null +++ b/src/extrapolate.F90 @@ -0,0 +1,62 @@ +! $Header$ + +! Extrapolate the difference between the grid function var and oldvar to +! var where indicated by mask. Use the normal direction given by dirx, +! diry, and dirz for extrapolation. + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "maskvalues.h" + +subroutine excision_extrapolate (ierr, cctkgh, var, oldvar, & + mask, dirx, diry, dirz, ni, nj, nk) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + ! Arguments + + ! out: zero for success, nonzero for error + integer :: ierr + + ! in: pointer to Cactus GH + CCTK_POINTER :: cctkgh + + ! in: array sizes for grid functions + ! (you can pass in cctk_lsh(:) for these) + integer :: ni,nj,nk + + ! inout: grid function that should be interpolated + CCTK_REAL :: var(ni,nj,nk) + + ! in: other grid function for interpolation + CCTK_REAL :: oldvar(ni,nj,nk) + + ! in: mask + CCTK_REAL :: mask(ni,nj,nk) + + ! in: normal directions to use for interpolation + CCTK_REAL :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk) + + integer i,j,k + integer ii,jj,kk + + do k=2,nk-1 + do j=2,nj-1 + do i=2,ni-1 + if (mask(i,j,k)==MASK_BOUNDARY) then + ii = i + dirx(i,j,k) + jj = j + diry(i,j,k) + kk = k + dirz(i,j,k) + var(i,j,k) = oldvar(i,j,k) + var(ii,jj,kk) - oldvar(ii,jj,kk) + end if + end do + end do + end do + + ierr = 0 + +end subroutine excision_extrapolate diff --git a/src/findboundary.F90 b/src/findboundary.F90 new file mode 100644 index 0000000..1e16e97 --- /dev/null +++ b/src/findboundary.F90 @@ -0,0 +1,57 @@ +! $Header$ + +! Take a mask that contains only the values 0 and 1. +! Return a mask where the outermost 0s have been replaced by 0.5s. + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "maskvalues.h" + +subroutine excision_findboundary (ierr, cctkgh, mask, ni, nj, nk) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + ! Arguments + + ! out: zero for success, nonzero for error + integer :: ierr + + ! in: pointer to Cactus GH + CCTK_POINTER :: cctkgh + + ! in: array sizes for grid functions + ! (you can pass in cctk_lsh(:) for these) + integer :: ni,nj,nk + + ! inout: mask + CCTK_REAL :: mask(ni,nj,nk) + + integer i,j,k + integer ii,jj,kk + logical bnd + + do k=2,nk-1 + do j=2,nj-1 + do i=2,ni-1 + if (mask(i,j,k)==MASK_EXCISED) then + bnd = .false. + do kk=k-1,k+1 + do jj=j-1,j+1 + do ii=i-1,i+1 + bnd = bnd .or. mask(i,j,k)/=MASK_EXCISED + end do + end do + end do + if (bnd) mask(i,j,k) = MASK_BOUNDARY + end if + end do + end do + end do + + ierr = 0 + +end subroutine excision_findboundary diff --git a/src/findnormals.F90 b/src/findnormals.F90 new file mode 100644 index 0000000..df0ae69 --- /dev/null +++ b/src/findnormals.F90 @@ -0,0 +1,290 @@ +! $Header$ + +! Take a mask where the boundary has been been marked. +! Return information about the normals in these locations. + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "maskvalues.h" + +subroutine excision_findnormals (ierr, cctkgh, mask, & + dirx, diry, dirz, ni, nj, nk) + + implicit none + + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + + ! Arguments + + ! out: zero for success, nonzero for error + integer :: ierr + + ! in: pointer to Cactus GH + CCTK_POINTER :: cctkgh + + ! in: array sizes for grid functions + ! (you can pass in cctk_lsh(:) for these) + integer :: ni,nj,nk + + ! in: mask + CCTK_REAL :: mask(ni,nj,nk) + + ! out: normal directions to use for interpolation + CCTK_REAL :: dirx(ni,nj,nk), diry(ni,nj,nk), dirz(ni,nj,nk) + + + + ! distinguish between faces, edges, and corners. + + ! look at the excised points, and distinguish the following cases, + ! in this order: + + ! [trivial case] + ! 1: if (none) + ! case 0: interior + ! [box shape cases] + ! 2: if (only one corner) + ! case I: corner + ! 3: if (only points on one edge) + ! case II: edge + ! 4: if (only points on one face) + ! case III: face + ! [convex cases] + ! 5: if (only points on two adjacent faces) + ! case IIa: anti-edge + ! 6: if (one points on three adjacent faces) + ! case Ia: anti-corner + ! [illegal cases] + ! error + + logical mask_interior (-1:+1,-1:+1,-1:+1, 1) + logical mask_corner (-1:+1,-1:+1,-1:+1, 8) + logical mask_edge (-1:+1,-1:+1,-1:+1, 12) + logical mask_face (-1:+1,-1:+1,-1:+1, 6) + logical mask_antiedge (-1:+1,-1:+1,-1:+1, 12) + logical mask_anticorner(-1:+1,-1:+1,-1:+1, 8) + integer dir_interior (3, 1), num_interior + integer dir_corner (3, 8), num_corner + integer dir_edge (3, 12), num_edge + integer dir_face (3, 6), num_face + integer dir_antiedge (3, 12), num_antiedge + integer dir_anticorner (3, 8), num_anticorner + + logical msk(-1:+1,-1:+1,-1:+1) + + integer di,dj,dk + integer i,j,k + integer n + integer dir(3) + + character(1000) message + + + + ! + ! Create templates + ! + num_interior = 0 + num_corner = 0 + num_edge = 0 + num_face = 0 + num_antiedge = 0 + num_anticorner = 0 + + do dk=-1,+1 + do dj=-1,+1 + do di=-1,+1 + + select case (count( (/di,dj,dk/) /= 0 )) + + case (0) + ! interior (or error, which is not explicitely checked) + num_interior = num_interior + 1 + dir_interior(:,num_interior) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_interior(i,j,k,num_interior) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case (1) + ! face (or anti-face, which is the same thing) + num_face = num_face + 1 + dir_face(:,num_face) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_face(i,j,k,num_face) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case (2) + ! edge or anti-edge + num_edge = num_edge + 1 + dir_edge(:,num_edge) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_edge(i,j,k,num_edge) = & + all( (/di,dj,dk/)==0 & + .or. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + num_antiedge = num_antiedge + 1 + dir_antiedge(:,num_antiedge) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_antiedge(i,j,k,num_antiedge) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case (3) + ! corner or anti-corner + num_corner = num_corner + 1 + dir_corner(:,num_corner) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_corner(i,j,k,num_corner) = & + all( (/di,dj,dk/)==0 & + .or. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + num_anticorner = num_anticorner + 1 + dir_anticorner(:,num_anticorner) = (/di,dj,dk/) + do k=-1,+1 + do j=-1,+1 + do i=-1,+1 + mask_anticorner(i,j,k,num_anticorner) = & + any( (/di,dj,dk/)/=0 & + .and. (/i,j,k/)==(/di,dj,dk/) ) + end do + end do + end do + + case default + call CCTK_WARN (0, "internal error") + end select + + end do + end do + end do + + ! sanity check + if ( num_interior /= 1 & + .or. num_corner /= 8 & + .or. num_edge /= 12 & + .or. num_face /= 6 & + .or. num_antiedge /= 12 & + .or. num_anticorner /= 8 ) then + call CCTK_WARN (0, "internal error") + end if + + + + ! + ! Calculate normal directions + ! + + dirx = 0 + diry = 0 + dirz = 0 + + do k=2,nk-1 + do j=2,nj-1 + do i=2,ni-1 + if (mask(i,j,k)==MASK_BOUNDARY) then + + ! make a snapshot + msk = mask(i-1:i+1,j-1:j+1,k-1:k+1) /= 0 + + ! for a template to match, all grid points have to be + ! either active, or set in the template. + search: do + + ! 1: case 0: interior + do n=1,num_interior + if (all(msk .or. mask_interior(:,:,:,n))) then + dir = dir_interior(:,n) + exit search + end if + end do + + ! 2: case I: corner + do n=1,num_corner + if (all(msk .or. mask_corner(:,:,:,n))) then + dir = dir_corner(:,n) + exit search + end if + end do + + ! 3: case II: edge + do n=1,num_edge + if (all(msk .or. mask_edge(:,:,:,n))) then + dir = dir_edge(:,n) + exit search + end if + end do + + ! 4: case III: face + do n=1,num_face + if (all(msk .or. mask_face(:,:,:,n))) then + dir = dir_face(:,n) + exit search + end if + end do + + ! 5: case IIa: anti-edge + do n=1,num_antiedge + if (all(msk .or. mask_antiedge(:,:,:,n))) then + dir = dir_antiedge(:,n) + exit search + end if + end do + + ! 6: case Ia: anticorner + do n=1,num_anticorner + if (all(msk .or. mask_anticorner(:,:,:,n))) then + dir = dir_anticorner(:,n) + exit search + end if + end do + + ! error: unclassified + call CCTK_WARN (0, "Mask is not convex") + + end do search + + dir(:) = -dir(:) + + dirx(i,j,k) = dir(1) + diry(i,j,k) = dir(2) + dirz(i,j,k) = dir(3) + + end if + + end do + end do + end do + + ! Success + ierr = 0 + +end subroutine excision_findnormals diff --git a/src/make.code.defn b/src/make.code.defn index af22b60..854a5e8 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = excision.F90 +SRCS = extrapolate.F90 findboundary.F90 findnormals.F90 # Subdirectories containing source files SUBDIRS = diff --git a/src/maskvalues.h b/src/maskvalues.h new file mode 100644 index 0000000..c2822f2 --- /dev/null +++ b/src/maskvalues.h @@ -0,0 +1,8 @@ +/* -*-f90-*- */ +/* $Header$ */ + +/* Values of the mask */ + +#define MASK_EXCISED 0 +#define MASK_BOUNDARY 0.5 +#define MASK_ACTIVE 1 -- cgit v1.2.3