Skip to content

Commit

Permalink
Merge pull request #266 from hannesbrandt/feature-more3d
Browse files Browse the repository at this point in the history 
Progress on 3d
  • Loading branch information
@donnaaboise
donnaaboise authored Jul 13, 2023
2 parents 90e0626 + fc17611 commit b79aae8
Show file tree
Hide file tree
Showing 4 changed files with 290 additions and 47 deletions.
129 changes: 91 additions & 38 deletions src/forestclaw2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -619,7 +619,7 @@ fclaw2d_patch_transform_face (fclaw2d_patch_t * ipatch,
#endif
)
{
double Rmx;
double iwidth;

FCLAW_ASSERT (ipatch->level == opatch->level);
FCLAW_ASSERT (0 <= ipatch->level && ipatch->level < P4EST_MAXLEVEL);
Expand All @@ -632,9 +632,9 @@ fclaw2d_patch_transform_face (fclaw2d_patch_t * ipatch,
FCLAW_ASSERT (opatch->zlower >= 0. && opatch->zlower < 1.);
#endif

FCLAW_ASSERT (mx >= 1 && mx == my);
FCLAW_ASSERT (mx >= 1 && my >= 1);
#ifdef P4_TO_P8
FCLAW_ASSERT (mx == mz);
FCLAW_ASSERT (mz >= 1);
#endif
FCLAW_ASSERT (based == 0 || based == 1);

Expand All @@ -650,15 +650,15 @@ fclaw2d_patch_transform_face (fclaw2d_patch_t * ipatch,
#endif

/* work with doubles -- exact for integers up to 52 bits of precision */
Rmx = (double) mx * (double) (1 << ipatch->level);
iwidth = (double) (1 << ipatch->level);

if (ftransform[8] & 4)
{
/* The two patches are in the same block. ftransform is not used */
*i += (int) ((ipatch->xlower - opatch->xlower) * Rmx);
*j += (int) ((ipatch->ylower - opatch->ylower) * Rmx);
*i += (int) ((ipatch->xlower - opatch->xlower) * iwidth * (double) mx);
*j += (int) ((ipatch->ylower - opatch->ylower) * iwidth * (double) my);
#ifdef P4_TO_P8
*k += (int) ((ipatch->zlower - opatch->zlower) * Rmx);
*k += (int) ((ipatch->zlower - opatch->zlower) * iwidth * (double) mz);
#endif
}
else
Expand All @@ -667,46 +667,72 @@ fclaw2d_patch_transform_face (fclaw2d_patch_t * ipatch,
const int *target_axis = &ftransform[3];
const int *edge_reverse = &ftransform[6];
double my_xyz[P4EST_DIM], target_xyz[P4EST_DIM];
double Rmxmymz[P4EST_DIM];
#ifdef FCLAW_ENABLE_DEBUG
int mxmymz[P4EST_DIM];

/* make mx, my and mz indexable */
mxmymz[0] = mx;
mxmymz[1] = my;
#ifdef P4_TO_P8
mxmymz[2] = mz;
#endif
#endif

/* make gridsize indexable */
Rmxmymz[0] = iwidth * (double) mx;
Rmxmymz[1] = iwidth * (double) my;
#ifdef P4_TO_P8
Rmxmymz[2] = iwidth * (double) mz;
#endif

/* the reference cube is stretched to mx times my units */
my_xyz[0] = ipatch->xlower * Rmx + *i - .5 * based;
my_xyz[1] = ipatch->ylower * Rmx + *j - .5 * based;
my_xyz[0] = ipatch->xlower * Rmxmymz[0] + *i - .5 * based;
my_xyz[1] = ipatch->ylower * Rmxmymz[1] + *j - .5 * based;
#ifdef P4_TO_P8
my_xyz[2] = ipatch->zlower * Rmx + *k - .5 * based;
my_xyz[2] = ipatch->zlower * Rmxmymz[2] + *k - .5 * based;
#endif

/* transform transversal directions */
FCLAW_ASSERT (mxmymz[my_axis[0]] == mxmymz[target_axis[0]]);
target_xyz[target_axis[0]] =
!edge_reverse[0] ? my_xyz[my_axis[0]] : Rmx - my_xyz[my_axis[0]];
!edge_reverse[0] ? my_xyz[my_axis[0]] : Rmxmymz[my_axis[0]] -
my_xyz[my_axis[0]];
#ifdef P4_TO_P8
FCLAW_ASSERT (mxmymz[my_axis[1]] == mxmymz[target_axis[1]]);
target_xyz[target_axis[1]] =
!edge_reverse[1] ? my_xyz[my_axis[1]] : Rmx - my_xyz[my_axis[1]];
!edge_reverse[1] ? my_xyz[my_axis[1]] : Rmxmymz[my_axis[1]] -
my_xyz[my_axis[1]];
#endif

/* transform normal direction */
FCLAW_ASSERT (mxmymz[my_axis[2]] == mxmymz[target_axis[2]]);
switch (edge_reverse[2])
{
case 0:
target_xyz[target_axis[2]] = -my_xyz[my_axis[2]];
break;
case 1:
target_xyz[target_axis[2]] = my_xyz[my_axis[2]] + Rmx;
target_xyz[target_axis[2]] =
my_xyz[my_axis[2]] + Rmxmymz[my_axis[2]];
break;
case 2:
target_xyz[target_axis[2]] = my_xyz[my_axis[2]] - Rmx;
target_xyz[target_axis[2]] =
my_xyz[my_axis[2]] - Rmxmymz[my_axis[2]];
break;
case 3:
target_xyz[target_axis[2]] = 2. * Rmx - my_xyz[my_axis[2]];
target_xyz[target_axis[2]] =
2. * Rmxmymz[my_axis[2]] - my_xyz[my_axis[2]];
break;
default:
SC_ABORT_NOT_REACHED ();
}

/* transform back to integer coordinates: this is exact */
*i = (int) (target_xyz[0] - opatch->xlower * Rmx + .5 * based);
*j = (int) (target_xyz[1] - opatch->ylower * Rmx + .5 * based);
*i = (int) (target_xyz[0] - opatch->xlower * Rmxmymz[0] + .5 * based);
*j = (int) (target_xyz[1] - opatch->ylower * Rmxmymz[1] + .5 * based);
#ifdef P4_TO_P8
*k = (int) (target_xyz[2] - opatch->zlower * Rmx + .5 * based);
*k = (int) (target_xyz[2] - opatch->zlower * Rmxmymz[2] + .5 * based);
#endif
}

Expand Down Expand Up @@ -736,7 +762,7 @@ fclaw2d_patch_transform_face2 (fclaw2d_patch_t * ipatch,
#ifdef P4_TO_P8
int dk;
#endif
double Rmx;
double owidth;

FCLAW_ASSERT (ipatch->level + 1 == opatch->level);
FCLAW_ASSERT (0 <= ipatch->level && opatch->level < P4EST_MAXLEVEL);
Expand All @@ -749,9 +775,9 @@ fclaw2d_patch_transform_face2 (fclaw2d_patch_t * ipatch,
FCLAW_ASSERT (opatch->zlower >= 0. && opatch->zlower < 1.);
#endif

FCLAW_ASSERT (mx >= 1 && mx == my);
FCLAW_ASSERT (mx >= 1 && my >= 1);
#ifdef P4_TO_P8
FCLAW_ASSERT (mx == mz);
FCLAW_ASSERT (mz >= 1);
#endif
FCLAW_ASSERT (based == 0 || based == 1);

Expand All @@ -767,20 +793,23 @@ fclaw2d_patch_transform_face2 (fclaw2d_patch_t * ipatch,
#endif

/* work with doubles -- exact for integers up to 52 bits of precision */
Rmx = (double) mx * (double) (1 << opatch->level);
owidth = (double) (1 << opatch->level);

if (ftransform[8] & 4)
{
int kx, ky, kz;

/* The two patches are in the same block. ftransform is undefined */
di = based + (int)
((ipatch->xlower - opatch->xlower) * Rmx + 2. * (*i - based));
dj = based + (int)
((ipatch->ylower - opatch->ylower) * Rmx + 2. * (*j - based));
((ipatch->xlower - opatch->xlower) * owidth * (double) mx +
2. * (*i - based));
dj = based +
(int) ((ipatch->ylower - opatch->ylower) * owidth * (double) my +
2. * (*j - based));
#ifdef P4_TO_P8
dk = based + (int)
((ipatch->zlower - opatch->zlower) * Rmx + 2. * (*k - based));
((ipatch->zlower - opatch->zlower) * owidth * (double) mz +
2. * (*k - based));
#else
kz = 0;
#endif
Expand Down Expand Up @@ -817,22 +846,43 @@ fclaw2d_patch_transform_face2 (fclaw2d_patch_t * ipatch,
const int *target_axis = &ftransform[3];
const int *edge_reverse = &ftransform[6];
double my_xyz[P4EST_DIM], target_xyz[P4EST_DIM];
double Rmxmymz[P4EST_DIM];
#ifdef FCLAW_ENABLE_DEBUG
int mxmymz[P4EST_DIM];

/* make mx, my and mz indexable */
mxmymz[0] = mx;
mxmymz[1] = my;
#ifdef P4_TO_P8
mxmymz[2] = mz;
#endif
#endif
/* make gridsize indexable */
Rmxmymz[0] = owidth * (double) mx;
Rmxmymz[1] = owidth * (double) my;
#ifdef P4_TO_P8
Rmxmymz[2] = owidth * (double) mz;
#endif

/* the reference cube is stretched to mx times my units */
my_xyz[0] = ipatch->xlower * Rmx + 2. * (*i + .5 - based);
my_xyz[1] = ipatch->ylower * Rmx + 2. * (*j + .5 - based);
my_xyz[0] = ipatch->xlower * Rmxmymz[0] + 2. * (*i + .5 - based);
my_xyz[1] = ipatch->ylower * Rmxmymz[1] + 2. * (*j + .5 - based);
#ifdef P4_TO_P8
my_xyz[2] = ipatch->zlower * Rmx + 2. * (*k + .5 - based);
my_xyz[2] = ipatch->zlower * Rmxmymz[2] + 2. * (*k + .5 - based);
#else
is[2] = ik[1] = ib[1] = in[2] = 0;
#endif

/* transform transversal directions */
FCLAW_ASSERT (mxmymz[target_axis[0]] == mxmymz[my_axis[0]]);
target_xyz[target_axis[0]] =
!edge_reverse[0] ? my_xyz[my_axis[0]] : Rmx - my_xyz[my_axis[0]];
!edge_reverse[0] ? my_xyz[my_axis[0]] : Rmxmymz[my_axis[0]] -
my_xyz[my_axis[0]];
#ifdef P4_TO_P8
FCLAW_ASSERT (mxmymz[target_axis[1]] == mxmymz[my_axis[1]]);
target_xyz[target_axis[1]] =
!edge_reverse[1] ? my_xyz[my_axis[1]] : Rmx - my_xyz[my_axis[1]];
!edge_reverse[1] ? my_xyz[my_axis[1]] : Rmxmymz[my_axis[1]] -
my_xyz[my_axis[1]];
#endif

/* transform normal direction */
Expand All @@ -842,23 +892,26 @@ fclaw2d_patch_transform_face2 (fclaw2d_patch_t * ipatch,
target_xyz[target_axis[2]] = -my_xyz[my_axis[2]];
break;
case 1:
target_xyz[target_axis[2]] = my_xyz[my_axis[2]] + Rmx;
target_xyz[target_axis[2]] =
my_xyz[my_axis[2]] + Rmxmymz[my_axis[2]];
break;
case 2:
target_xyz[target_axis[2]] = my_xyz[my_axis[2]] - Rmx;
target_xyz[target_axis[2]] =
my_xyz[my_axis[2]] - Rmxmymz[my_axis[2]];
break;
case 3:
target_xyz[target_axis[2]] = 2. * Rmx - my_xyz[my_axis[2]];
target_xyz[target_axis[2]] =
2. * Rmxmymz[my_axis[2]] - my_xyz[my_axis[2]];
break;
default:
SC_ABORT_NOT_REACHED ();
}

/* move back into integer coordinates: this is exact */
di = (int) (target_xyz[0] - opatch->xlower * Rmx) + based - 1;
dj = (int) (target_xyz[1] - opatch->ylower * Rmx) + based - 1;
di = (int) (target_xyz[0] - opatch->xlower * Rmxmymz[0]) + based - 1;
dj = (int) (target_xyz[1] - opatch->ylower * Rmxmymz[1]) + based - 1;
#ifdef P4_TO_P8
dk = (int) (target_xyz[2] - opatch->zlower * Rmx) + based - 1;
dk = (int) (target_xyz[2] - opatch->zlower * Rmxmymz[2]) + based - 1;
#endif

/* Run through the child cells in order of the small patch */
Expand Down Expand Up @@ -969,7 +1022,7 @@ fclaw2d_patch_corner_neighbors (fclaw2d_domain_t * domain,
FCLAW_ASSERT (0 <= qid);
if (qid >= mesh->local_num_quadrants + mesh->ghost_num_quadrants)
{
/* This is an inter-tree (face or corner) corner neighbor */
/* This is an inter-tree (face, edge or corner) corner neighbor */
cornerid =
qid - (mesh->local_num_quadrants + mesh->ghost_num_quadrants);
FCLAW_ASSERT (cornerid < mesh->local_num_corners);
Expand Down
6 changes: 4 additions & 2 deletions src/forestclaw2d.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -507,7 +507,8 @@ int fclaw2d_patch_face_transformation_valid (const int ftransform[]);
* we require \a ftransform[8] |= 4.
* \param [in] mx Number of cells along x direction of patch.
* \param [in] my Number of cells along y direction of patch.
* This function assumes \a mx == \a my.
* The number of cells must match according to the face
* transformation.
* \param [in] based Indices are 0-based for corners and 1-based for cells.
* \param [in,out] i Integer coordinate along x-axis in \a based .. \a mx.
* \param [in,out] j Integer coordinate along y-axis in \a based .. \a my.
Expand All @@ -530,7 +531,8 @@ void fclaw2d_patch_transform_face (fclaw2d_patch_t * ipatch,
* we require \a ftransform[8] |= 4.
* \param [in] mx Number of cells along x direction of patch.
* \param [in] my Number of cells along y direction of patch.
* This function assumes \a mx == \a my.
* The number of cells must match according to the face
* transformation.
* \param [in] based Indices are 0-based for corners and 1-based for cells.
* \param [in,out] i FOUR (4) integer coordinates along x-axis in
* \a based .. \a mx. On input, only the first is used.
Expand Down
Loading

0 comments on commit b79aae8

Please sign in to comment.