From 6883a7e95eca9c8eeb3b3833ceb3404cd536894f Mon Sep 17 00:00:00 2001 From: schnetter Date: Thu, 4 Jan 2007 15:40:06 +0000 Subject: Add support for 4-tensors. Patch mostly thanks to Bela Szilagyi. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@19 082bdb00-0f4f-0410-b49e-b1835e5f2039 --- src/apply.c | 170 +++++++++++---------- src/interpolate.c | 431 ++++++++++++++++++++++++++++-------------------------- 2 files changed, 304 insertions(+), 297 deletions(-) diff --git a/src/apply.c b/src/apply.c index 1f71e1b..918957e 100644 --- a/src/apply.c +++ b/src/apply.c @@ -47,12 +47,9 @@ CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_apply_c); assert(jjmax>=-1&&jjmin=-1&&kkmin= CCTK_NumVars()) - { + if (vi < 0 || vi >= CCTK_NumVars()) { CCTK_WARN (0, "Illegal variable index"); } - if (verbose) - { + if (verbose) { fullname = CCTK_FullName (vi); - if (! fullname) - { + if (! fullname) { CCTK_WARN (0, "Internal error in CCTK_FullName"); } CCTK_VInfo (CCTK_THORNSTRING, @@ -197,8 +189,7 @@ BndReflectVI (cGH const * restrict const cctkGH, /* Get and check group information */ gi = CCTK_GroupIndexFromVarI (vi); - if (gi < 0 || gi > CCTK_NumGroups()) - { + if (gi < 0 || gi > CCTK_NumGroups()) { CCTK_WARN (0, "Internal error in CCTK_GroupIndexFromVarI"); } @@ -258,11 +249,14 @@ BndReflectVI (cGH const * restrict const cctkGH, ttype = UNKNOWN; tcomponent = 0; - if (CCTK_EQUALS (tensortypealias, "scalar")) - { + if (CCTK_EQUALS (tensortypealias, "scalar")) { /* scalar */ ttype = SCALAR; tcomponent = 0; + } else if (CCTK_EQUALS (tensortypealias, "4scalar")) { + /* 4-scalar */ + ttype = SCALAR; + tcomponent = 0; } else if (CCTK_EQUALS (tensortypealias, "u") || CCTK_EQUALS (tensortypealias, "d")) { @@ -270,12 +264,40 @@ BndReflectVI (cGH const * restrict const cctkGH, assert (numvars == 3); ttype = VECTOR; tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "4u") + || CCTK_EQUALS (tensortypealias, "4d")) + { + /* 4-vector */ + assert (numvars == 4); + if (vi == firstvar) { + ttype = SCALAR; + tcomponent = 0; + } else { + ttype = VECTOR; + tcomponent = vi - firstvar - 1; + } } else if (CCTK_EQUALS (tensortypealias, "uu_sym") - || CCTK_EQUALS (tensortypealias, "dd_sym")) { + || CCTK_EQUALS (tensortypealias, "dd_sym")) + { /* symmetric tensor */ assert (numvars == 6); ttype = TENSOR; tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "4uu_sym") + || CCTK_EQUALS (tensortypealias, "4dd_sym")) + { + /* symmetric 4-tensor */ + assert (numvars == 10); + if (vi == firstvar) { + ttype = SCALAR; + tcomponent = 0; + } else if (vi <= firstvar+3) { + ttype = VECTOR; + tcomponent = vi - firstvar - 1; + } else { + ttype = TENSOR; + tcomponent = vi - firstvar - 4; + } } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) { /* Weyl scalars, stored as 10 real values. NOTE: This assumes that Psi_0 comes first, which is NOT the default with @@ -292,8 +314,7 @@ BndReflectVI (cGH const * restrict const cctkGH, free (groupname); } - switch (ttype) - { + switch (ttype) { case SCALAR: assert (tcomponent>=0 && tcomponent<1); break; @@ -342,21 +363,16 @@ BndReflectVI (cGH const * restrict const cctkGH, /* Loop over all directions and faces */ - for (dir=0; dir<3; ++dir) - { - for (face=0; face<2; ++face) - { + for (dir=0; dir<3; ++dir) { + for (face=0; face<2; ++face) { /* If there is a reflection symmetry on that face */ - if (do_reflection[2*dir+face]) - { + if (do_reflection[2*dir+face]) { /* If we have the outer boundary of that face */ - if (cctkGH->cctk_bbox[2*dir+face]) - { + if (cctkGH->cctk_bbox[2*dir+face]) { /* Find parity */ parity = tensorparity; - switch (ttype) - { + switch (ttype) { case SCALAR: parity *= +1; break; @@ -364,8 +380,7 @@ BndReflectVI (cGH const * restrict const cctkGH, parity *= dir == tcomponent ? -1 : +1; break; case TENSOR: - switch (tcomponent) - { + switch (tcomponent) { case 0: parity *= +1; break; case 1: parity *= (dir == 0 || dir == 1) ? -1 : +1; break; case 2: parity *= (dir == 0 || dir == 2) ? -1 : +1; break; @@ -375,44 +390,39 @@ BndReflectVI (cGH const * restrict const cctkGH, default: assert (0); } break; - case WEYLSCALARS_REAL: - { - static int const weylparities[10][3] = - {{+1,+1,+1}, - {-1,-1,-1}, - {+1,+1,+1}, - {-1,-1,-1}, - {+1,+1,+1}, - {-1,-1,-1}, - {+1,+1,+1}, - {-1,-1,-1}, + case WEYLSCALARS_REAL: { + static int const weylparities[10][3] = + {{+1,+1,+1}, + {-1,-1,-1}, + {+1,+1,+1}, + {-1,-1,-1}, + {+1,+1,+1}, + {-1,-1,-1}, + {+1,+1,+1}, + {-1,-1,-1}, {+1,+1,+1}, - {-1,-1,-1}}; - parity *= weylparities[tcomponent][dir]; - } + {-1,-1,-1}}; + parity *= weylparities[tcomponent][dir]; break; + } default: assert (0); } /* Find region extent */ - for (d=0; d<3; ++d) - { + for (d=0; d<3; ++d) { lsh[d] = cctkGH->cctk_lsh[d]; imin[d] = 0; imax[d] = cctkGH->cctk_lsh[d]; ioff[d] = 0; idir[d] = 1; } - if (face == 0) - { + if (face == 0) { imax[dir] = cctkGH->cctk_nghostzones[dir]; ioff[dir] = (+ 2*cctkGH->cctk_nghostzones[dir] - 1 + (do_stagger[2*dir+face] ? 0 : 1)); idir[dir] = -1; - } - else - { + } else { imin[dir] = cctkGH->cctk_lsh[dir] - cctkGH->cctk_nghostzones[dir]; ioff[dir] = (- 2*cctkGH->cctk_nghostzones[dir] - (do_stagger[2*dir+face] ? 0 : 1)); @@ -420,8 +430,7 @@ BndReflectVI (cGH const * restrict const cctkGH, } /* Copy region */ - switch (group.vartype) - { + switch (group.vartype) { #ifdef HAVE_CCTK_INT1 case CCTK_VARIABLE_INT1: @@ -725,86 +734,71 @@ ReflectionSymmetry_Apply (CCTK_ARGUMENTS) int istat; int ierr; - if (! cctkGH) - { + if (! cctkGH) { CCTK_WARN (0, "Argument cctkGH is NULL"); } nvars = Boundary_SelectedGVs (cctkGH, 0, NULL, NULL, NULL, NULL, NULL); - if (nvars < 0) - { + if (nvars < 0) { CCTK_WARN (0, "Internal error in Boundary_SelectedGVs"); } - if (nvars == 0) - { + if (nvars == 0) { /* Nothing to do */ return; } indices = malloc (nvars * sizeof *indices); - if (! indices) - { + if (! indices) { CCTK_WARN (0, "Out of memory"); } faces = malloc (nvars * sizeof *faces); - if (! faces) - { + if (! faces) { CCTK_WARN (0, "Out of memory"); } widths = malloc (nvars * sizeof *widths); - if (! widths) - { + if (! widths) { CCTK_WARN (0, "Out of memory"); } tables = malloc (nvars * sizeof *tables); - if (! tables) - { + if (! tables) { CCTK_WARN (0, "Out of memory"); } istat = Boundary_SelectedGVs (cctkGH, nvars, indices, faces, widths, tables, 0); - if (istat != nvars) - { + if (istat != nvars) { CCTK_WARN (0, "Internal error in Boundary_SelectedGVs"); } - for (i=0; i= CCTK_NumVars()) - { + if (vi < 0 || vi >= CCTK_NumVars()) { CCTK_WARN (0, "Illegal variable index"); } - if (widths[i] < 0) - { + if (widths[i] < 0) { CCTK_WARN (0, "Illegal boundary width"); } dim = CCTK_GroupDimFromVarI (vi); - if (dim < 0) - { + if (dim < 0) { CCTK_WARN (0, "Illegal dimension"); } stencil = malloc (dim * sizeof *stencil); - if (! stencil) - { + if (! stencil) { CCTK_WARN (0, "Out of memory"); } ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi); - if (ierr) - { + if (ierr) { CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI"); } CheckBoundaryParameters (cctkGH, vi, stencil); ierr = BndReflectVI (cctkGH, stencil, vi); - if (ierr) - { + if (ierr) { CCTK_WARN (0, "Internal error in BndReflectVI"); } diff --git a/src/interpolate.c b/src/interpolate.c index e4c242d..eb60460 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -67,10 +67,8 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, do_reflection[5] = reflection_upper_z; newfaces = faces; - for (d=0; d<6; ++d) - { - if (do_reflection[d]) - { + for (d=0; d<6; ++d) { + if (do_reflection[d]) { assert (newfaces & (1 << d)); newfaces &= ~ (1 << d); } @@ -80,25 +78,20 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, /* Fold coordinates */ assert (interp_coords_type == CCTK_VARIABLE_REAL); - for (dir=0; dir<3; ++dir) - { + for (dir=0; dir<3; ++dir) { assert (! do_reflection[2*dir+1]); - if (do_reflection[2*dir]) - { + if (do_reflection[2*dir]) { new_interp_coords[dir] = malloc (N_interp_points * sizeof (CCTK_REAL)); assert (N_interp_points == 0 || new_interp_coords[dir]); - for (n=0; n=0 && vi=0 && gi0); - firstvar = CCTK_FirstVarIndexI(gi); - assert (firstvar>=0); - table = CCTK_GroupTagsTableI(gi); - assert (table>=0); - - - - /* Get and check tensor type information */ - ierr = Util_TableGetString - (table, sizeof tensortypealias, tensortypealias, "tensortypealias"); - if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - groupname = CCTK_GroupName(gi); - assert (groupname); - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Tensor type alias not declared for group \"%s\" -- assuming a scalar", - groupname); - free (groupname); - strcpy (tensortypealias, "scalar"); - } else if (ierr<0) { - groupname = CCTK_GroupName(gi); - assert (groupname); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Error in tensor type alias declaration for group \"%s\"", - groupname); - free (groupname); - } - - ttype = UNKNOWN; - tcomponent = 0; - if (CCTK_EQUALS (tensortypealias, "scalar")) { - /* scalar */ - if (numvars != 1) { + for (m=0; m=0 && vi=0 && gi0); + firstvar = CCTK_FirstVarIndexI(gi); + assert (firstvar>=0); + table = CCTK_GroupTagsTableI(gi); + assert (table>=0); + + + + /* Get and check tensor type information */ + ierr = Util_TableGetString + (table, sizeof tensortypealias, tensortypealias, "tensortypealias"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Tensor type alias not declared for group \"%s\" -- assuming a scalar", + groupname); + free (groupname); + strcpy (tensortypealias, "scalar"); + } else if (ierr<0) { groupname = CCTK_GroupName(gi); assert (groupname); - CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element", + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in tensor type alias declaration for group \"%s\"", groupname); free (groupname); } - ttype = SCALAR; + + ttype = UNKNOWN; tcomponent = 0; - } - else if (CCTK_EQUALS (tensortypealias, "u") - || CCTK_EQUALS (tensortypealias, "d")) { - /* vector */ - assert (numvars == 3); - ttype = VECTOR; - tcomponent = vi - firstvar; - } else if (CCTK_EQUALS (tensortypealias, "uu_sym") - || CCTK_EQUALS (tensortypealias, "dd_sym")) { - /* symmetric tensor */ - assert (numvars == 6); - ttype = TENSOR; - tcomponent = vi - firstvar; - } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) { - /* Weyl scalars, stored as 10 real numbers */ - assert (numvars == 10); - ttype = WEYLSCALARS_REAL; - tcomponent = vi - firstvar; - } else { - char * groupname = CCTK_GroupName(gi); - assert (groupname); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Illegal tensor type alias for group \"%s\"", - groupname); - free (groupname); - } - - switch (ttype) - { - case SCALAR: - assert (tcomponent>=0 && tcomponent<1); - break; - case VECTOR: - assert (tcomponent>=0 && tcomponent<3); - break; - case TENSOR: - assert (tcomponent>=0 && tcomponent<6); - break; - case WEYLSCALARS_REAL: - assert (tcomponent>=0 && tcomponent<10); - break; - default: - assert (0); - } - - ierr = Util_TableGetInt (table, & tensorparity, "tensorparity"); - if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - tensorparity = +1; - } else if (ierr<0) { - groupname = CCTK_GroupName(gi); - assert (groupname); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Error in tensor parity declaration for group \"%s\"", - groupname); - free (groupname); - } - - - - /* Calculate parities */ - parities[0] = parities[1] = parities[2] = +1; - switch (ttype) - { - case SCALAR: - /* do nothing */ - break; - case VECTOR: - parities[tcomponent] = -1; - break; - case TENSOR: - switch (tcomponent) + if (CCTK_EQUALS (tensortypealias, "scalar")) { + /* scalar */ + if (numvars != 1) { + groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element", + groupname); + free (groupname); + } + ttype = SCALAR; + tcomponent = 0; + } else if (CCTK_EQUALS (tensortypealias, "4scalar")) { + /* 4-scalar */ + assert (numvars == 1); + ttype = SCALAR; + tcomponent = 0; + } else if (CCTK_EQUALS (tensortypealias, "u") + || CCTK_EQUALS (tensortypealias, "d")) { - case 0: break; - case 1: parities[0] = parities[1] = -1; break; - case 2: parities[0] = parities[2] = -1; break; - case 3: break; - case 4: parities[1] = parities[2] = -1; break; - case 5: break; - default: assert (0); - } - break; - case WEYLSCALARS_REAL: + /* vector */ + assert (numvars == 3); + ttype = VECTOR; + tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "4u") + || CCTK_EQUALS (tensortypealias, "4d")) + { + /* 4-vector */ + assert (numvars == 4); + if (vi == firstvar) { + ttype = SCALAR; + tcomponent = 0; + } else { + ttype = VECTOR; + tcomponent = vi - firstvar - 1; + } + } else if (CCTK_EQUALS (tensortypealias, "uu_sym") + || CCTK_EQUALS (tensortypealias, "dd_sym")) + { + /* symmetric tensor */ + assert (numvars == 6); + ttype = TENSOR; + tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "4uu_sym") + || CCTK_EQUALS (tensortypealias, "4dd_sym")) { + /* symmetric 4-tensor */ + assert (numvars == 10); + if (vi == firstvar) { + ttype = SCALAR; + tcomponent = 0; + } else if (vi <= firstvar+3) { + ttype = VECTOR; + tcomponent = vi - firstvar - 1; + } else { + ttype = TENSOR; + tcomponent = vi - firstvar - 4; + } + } else if (CCTK_EQUALS (tensortypealias, "weylscalars_real")) { + /* Weyl scalars, stored as 10 real numbers */ + assert (numvars == 10); + ttype = WEYLSCALARS_REAL; + tcomponent = vi - firstvar; + } else { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Illegal tensor type alias for group \"%s\"", + groupname); + free (groupname); + } + + switch (ttype) { + case SCALAR: + assert (tcomponent>=0 && tcomponent<1); + break; + case VECTOR: + assert (tcomponent>=0 && tcomponent<3); + break; + case TENSOR: + assert (tcomponent>=0 && tcomponent<6); + break; + case WEYLSCALARS_REAL: + assert (tcomponent>=0 && tcomponent<10); + break; + default: + assert (0); + } + + ierr = Util_TableGetInt (table, & tensorparity, "tensorparity"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + tensorparity = +1; + } else if (ierr<0) { + groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in tensor parity declaration for group \"%s\"", + groupname); + free (groupname); + } + + + + /* Calculate parities */ + parities[0] = parities[1] = parities[2] = +1; + switch (ttype) { + case SCALAR: + /* do nothing */ + break; + case VECTOR: + parities[tcomponent] = -1; + break; + case TENSOR: + switch (tcomponent) { + case 0: break; + case 1: parities[0] = parities[1] = -1; break; + case 2: parities[0] = parities[2] = -1; break; + case 3: break; + case 4: parities[1] = parities[2] = -1; break; + case 5: break; + default: assert (0); + } + break; + case WEYLSCALARS_REAL: { static int const weylparities[10][3] = {{+1,+1,+1}, {-1,-1,-1}, @@ -333,62 +352,56 @@ ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, for (dir=0; dir<3; ++dir) { parities[dir] = weylparities[tcomponent][dir]; } + break; } - break; - default: - assert (0); - } - - - - /* Take derivatives into account */ - { - int code = operation_codes[m]; - while (code) { - const int d = code % 10 - 1; - code /= 10; - assert (d>=0 && d<3); - parities[d] *= -1; + default: + assert (0); } - } - - - - /* Are there negative parities? */ - needs_checking = 0; - for (dir=0; dir<3; ++dir) - { - check_dir[dir] = do_reflection[2*dir] && parities[dir] < 0; - needs_checking |= check_dir[dir]; - } - - - - /* Loop over all points and unfold */ - if (needs_checking) - { - for (n=0; n=0 && d<3); + parities[d] *= -1; + } + } + + + + /* Are there negative parities? */ + needs_checking = 0; + for (dir=0; dir<3; ++dir) { + check_dir[dir] = do_reflection[2*dir] && parities[dir] < 0; + needs_checking |= check_dir[dir]; + } + + + + /* Loop over all points and unfold */ + if (needs_checking) { + for (n=0; n