aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetReduce/src/mask_coords.c
blob: d5d8df7c4062cadabb760bb70da9536096154289 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/mask_coords.c,v 1.3 2004/08/04 13:03:09 schnetter Exp $ */

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"



void
CoordBase_SetupMask (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  
  CCTK_INT nboundaryzones[6];
  CCTK_INT is_internal[6];
  CCTK_INT is_staggered[6];
  CCTK_INT shiftout[6];
  
  int imin[3], imax[3];         /* boundary extent */
  
  int i, j, k;
  int d, f;
  int dd;
  
  int ierr;
  
  
  
  ierr = GetBoundarySpecification
    (6, nboundaryzones, is_internal, is_staggered, shiftout);
  if (ierr != 0) {
    CCTK_WARN (0, "Could not get boundary specification");
  }
  
  /* Loop over all dimensions and faces */
  for (d=0; d<3; ++d) {
    for (f=0; f<2; ++f) {
      /* If this processor has the outer boundary */
      if (cctk_bbox[2*d+f]) {
        
        int npoints;
        
        if (is_internal[2*d+f]) {
          /* The boundary extends inwards */
          npoints = - nboundaryzones[2*d+f] + shiftout[2*d+f];
        } else {
          /* The boundary extends outwards */
          npoints = nboundaryzones[2*d+f] + shiftout[2*d+f] - 1;
        }
        
        /* If there are boundary that should be ignored */
        if (npoints >= 0) {
          
          if (npoints < 0 || npoints > cctk_lsh[d]) {
            CCTK_WARN (0, "Illegal number of boundary points");
          }
          
          /* Calculate the extent of the boundary */
          for (dd=0; dd<3; ++dd) {
            imin[dd] = 0;
            imax[dd] = cctk_lsh[dd];
          }
          
          if (f==0) {
            /* lower face */
            imax[d] = imin[d] + npoints;
          } else {
            /* upper face */
            imin[d] = imax[d] - npoints;
          }
          
          /* Loop over the boundary */
          if (verbose) {
            CCTK_VInfo (CCTK_THORNSTRING,
                        "Setting boundary points in direction %d face %d to weight 0", d, f);
          }
          for (k=imin[2]; k<imax[2]; ++k) {
            for (j=imin[1]; j<imax[1]; ++j) {
              for (i=imin[0]; i<imax[0]; ++i) {
                
                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
                weight[ind] = 0.0;
                
              }
            }
          }
          
          /* When the boundary is not staggered, then give the points
             on the boundary the weight 1/2 */
          if (! is_staggered[2*d+f]) {
            
            /* Adapt the extent of the boundary */
            if (f==0) {
              /* lower face */
              imin[d] = imax[d];
              imax[d] = imin[d] + 1;
            } else {
              /* upper face */
              imax[d] = imin[d];
              imin[d] = imax[d] - 1;
            }
          
            /* Loop over the points next to boundary */
            if (verbose) {
              CCTK_VInfo (CCTK_THORNSTRING,
                          "Setting staggered boundary points in direction %d face %d to weight 1/2", d, f);
            }
            for (k=imin[2]; k<imax[2]; ++k) {
              for (j=imin[1]; j<imax[1]; ++j) {
                for (i=imin[0]; i<imax[0]; ++i) {
                  
                  int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
                  weight[ind] *= 0.5;
                  
                }
              }
            }
            
          } /* if the boundary is not staggered */
          
        } /* if there are boundary points */
        
      } /* if is outer boundary */
    } /* loop over faces */
  } /* loop over directions */
  
}