From d02895af18ae9c3673202d6231da078cb3d74bac Mon Sep 17 00:00:00 2001 From: schnetter Date: Tue, 1 Jun 2004 16:09:35 +0000 Subject: A thorn that applies reflection symmetry boundary conditions, including symmetry interpolation. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@2 082bdb00-0f4f-0410-b49e-b1835e5f2039 --- README | 10 + doc/documentation.tex | 144 ++++++++++++ interface.ccl | 68 ++++++ param.ccl | 62 +++++ schedule.ccl | 12 + src/apply.c | 627 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/interpolate.c | 359 +++++++++++++++++++++++++++++ src/make.code.defn | 8 + src/reflection.h | 29 +++ src/register.c | 66 ++++++ 10 files changed, 1385 insertions(+) create mode 100644 README create mode 100644 doc/documentation.tex create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/apply.c create mode 100644 src/interpolate.c create mode 100644 src/make.code.defn create mode 100644 src/reflection.h create mode 100644 src/register.c diff --git a/README b/README new file mode 100644 index 0000000..0f56f7a --- /dev/null +++ b/README @@ -0,0 +1,10 @@ +CVS info : $Header$ + +Cactus Code Thorn ReflectionSymmetry +Thorn Author(s) : Erik Schnetter +Thorn Maintainer(s) : Erik Schnetter +-------------------------------------------------------------------------- + +Purpose of the thorn: + +Provide reflection symmetries, i.e. bitant, quadrant, and octant mode. diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..fbc1f83 --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,144 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header$ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevant thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date, etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be separated with a \\ or a comma. +% - You can define your own macros, but they must appear after +% the START CACTUS THORNGUIDE line, and must not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphicx package. +% More specifically, with the "\includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us to create a PDF version of the ThornGuide +% via pdflatex. +% - References should be included with the latex "\bibitem" command. +% - Use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - Do not use \appendix, instead include any appendices you need as +% standard sections. +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header$ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../../doc/latex/cactus} + +\begin{document} + +% The author of the documentation +\author{Erik Schnetter \textless schnetter@aei.mpg.de\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{} + +% the date your document was last changed, if your document is in CVS, +% please use: +% \date{$ $Date$ $} +\date{June 01 2004} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc + +% Add an abstract for this thorn's documentation +\begin{abstract} + +\end{abstract} + +% The following sections are suggestive only. +% Remove them or add your own. + +\section{Introduction} + +\section{Physical System} + +\section{Numerical Implementation} + +\section{Using This Thorn} + +\subsection{Obtaining This Thorn} + +\subsection{Basic Usage} + +\subsection{Special Behaviour} + +\subsection{Interaction With Other Thorns} + +\subsection{Examples} + +\subsection{Support and Feedback} + +\section{History} + +\subsection{Thorn Source Code} + +\subsection{Thorn Documentation} + +\subsection{Acknowledgements} + + +\begin{thebibliography}{9} + +\end{thebibliography} + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..1a335bc --- /dev/null +++ b/interface.ccl @@ -0,0 +1,68 @@ +# Interface definition for thorn ReflectionSymmetry +# $Header$ + +IMPLEMENTS: ReflectionSymmetry + + + +CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name) +REQUIRES FUNCTION SymmetryRegister + +CCTK_INT FUNCTION \ + SymmetryRegisterGrid \ + (CCTK_POINTER IN cctkGH, \ + CCTK_INT IN sym_handle, \ + CCTK_INT IN ARRAY which_faces, \ # array [N_FACES] + CCTK_INT IN ARRAY symmetry_zone_width) # array [N_FACES] +REQUIRES FUNCTION SymmetryRegisterGrid + +CCTK_INT FUNCTION \ + SymmetryRegisterGridInterpolator \ + (CCTK_POINTER IN cctkGH, \ + CCTK_INT IN sym_handle, \ + CCTK_INT CCTK_FPOINTER IN symmetry_interpolate \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN N_dims, \ + CCTK_INT IN local_interp_handle, \ + CCTK_INT IN param_table_handle, \ + CCTK_INT IN coord_system_handle, \ + CCTK_INT IN N_interp_points, \ + CCTK_INT IN interp_coords_type, \ + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \ + CCTK_INT IN N_input_arrays, \ + CCTK_INT ARRAY IN input_array_indices, \ + CCTK_INT IN N_output_arrays, \ + CCTK_INT ARRAY IN output_array_types, \ + CCTK_POINTER ARRAY IN output_arrays, \ + CCTK_INT IN faces)) +REQUIRES FUNCTION SymmetryRegisterGridInterpolator + +CCTK_INT FUNCTION \ + SymmetryInterpolateFaces \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN N_dims, \ + CCTK_INT IN local_interp_handle, \ + CCTK_INT IN param_table_handle, \ + CCTK_INT IN coord_system_handle, \ + CCTK_INT IN N_interp_points, \ + CCTK_INT IN interp_coords_type, \ + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \ + CCTK_INT IN N_input_arrays, \ + CCTK_INT ARRAY IN input_array_indices, \ + CCTK_INT IN N_output_arrays, \ + CCTK_INT ARRAY IN output_array_types, \ + CCTK_POINTER ARRAY IN output_arrays, \ + CCTK_INT IN faces) +REQUIRES FUNCTION SymmetryInterpolateFaces + + + +CCTK_INT FUNCTION Boundary_SelectedGVs \ + (CCTK_POINTER_TO_CONST IN GH, \ + CCTK_INT IN array_size, \ + CCTK_INT ARRAY OUT var_indicies, \ + CCTK_INT ARRAY OUT faces, \ + CCTK_INT ARRAY OUT boundary_widths, \ + CCTK_INT ARRAY OUT table_handles, \ + CCTK_STRING IN bc_name) +REQUIRES FUNCTION Boundary_SelectedGVs diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..f6406ec --- /dev/null +++ b/param.ccl @@ -0,0 +1,62 @@ +# Parameter definitions for thorn ReflectionSymmetry +# $Header$ + +BOOLEAN verbose "Produce screen output while applying boundary conditions" +{ +} "no" + + + +BOOLEAN reflection_x "Reflection symmetry at the lower x boundary" +{ +} "no" + +BOOLEAN reflection_y "Reflection symmetry at the lower y boundary" +{ +} "no" + +BOOLEAN reflection_z "Reflection symmetry at the lower z boundary" +{ +} "no" + + + +BOOLEAN reflection_upper_x "Reflection symmetry at the upper x boundary" +{ +} "no" + +BOOLEAN reflection_upper_y "Reflection symmetry at the upper y boundary" +{ +} "no" + +BOOLEAN reflection_upper_z "Reflection symmetry at the upper z boundary" +{ +} "no" + + + +BOOLEAN avoid_origin_x "Stagger about the origin on the lower x boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_y "Stagger about the origin on the lower y boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_z "Stagger about the origin on the lower z boundary?" +{ +} "yes" + + + +BOOLEAN avoid_origin_upper_x "Stagger about the origin on the upper x boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_upper_y "Stagger about the origin on the upper y boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_upper_z "Stagger about the origin on the upper z boundary?" +{ +} "yes" diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..8c5db0e --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,12 @@ +# Schedule definitions for thorn ReflectionSymmetry +# $Header$ + +SCHEDULE ReflectionSymmetry_Register IN SymmetryRegister +{ + LANG: C +} "Register reflection symmetry boundaries" + +SCHEDULE ReflectionSymmetry_Apply IN BoundaryConditions +{ + LANG: C +} "Apply reflection symmetries" diff --git a/src/apply.c b/src/apply.c new file mode 100644 index 0000000..114c1f7 --- /dev/null +++ b/src/apply.c @@ -0,0 +1,627 @@ +#include +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "reflection.h" + + + +static const char rcsid[] = "$Header$"; +CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_apply_c); + + + +#define COPY(VARTYPE) \ + static void \ + copy_##VARTYPE (VARTYPE const * restrict const srcvar, \ + VARTYPE * restrict const dstvar, \ + int const ni, int const nj, int const nk, \ + int const imin, int const jmin, int const kmin, \ + int const imax, int const jmax, int const kmax, \ + int const ioff, int const joff, int const koff, \ + int const idir, int const jdir, int const kdir, \ + int const parity) \ + { \ + int i, j, k; \ + \ + for (k=kmin; k= CCTK_NumVars()) + { + CCTK_WARN (0, "Illegal variable index"); + } + + if (verbose) + { + fullname = CCTK_FullName (vi); + if (! fullname) + { + CCTK_WARN (0, "Internal error in CCTK_FullName"); + } + CCTK_VInfo (CCTK_THORNSTRING, + "Applying reflection boundary conditions to \"%s\"", + fullname); + free (fullname); + } + + + + /* Get and check group information */ + gi = CCTK_GroupIndexFromVarI (vi); + if (gi < 0 || gi > CCTK_NumGroups()) + { + CCTK_WARN (0, "Internal error in CCTK_GroupIndexFromVarI"); + } + + ierr = CCTK_GroupData (gi, &group); + assert (!ierr); + assert (group.grouptype == CCTK_GF); + assert (group.disttype == CCTK_DISTRIB_DEFAULT); + assert (group.stagtype == 0); + + firstvar = CCTK_FirstVarIndexI (gi); + assert (firstvar>=0 && firstvar=0); + + ierr = CCTK_GroupDynamicData (cctkGH, gi, &data); + assert (!ierr); + + table = CCTK_GroupTagsTableI(gi); + assert (table>=0); + + varptr = CCTK_VarDataPtrI (cctkGH, 0, vi); + assert (varptr); + + + + /* 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) { + 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, "u")) + { + /* vector */ + assert (numvars == 3); + ttype = VECTOR; + tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) { + /* symmetric tensor */ + assert (numvars == 6); + ttype = TENSOR; + 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; + default: + assert (0); + } + + + + /* Reflection symmetry information */ + do_reflection[0] = reflection_x; + do_reflection[1] = reflection_upper_x; + do_reflection[2] = reflection_y; + do_reflection[3] = reflection_upper_y; + do_reflection[4] = reflection_z; + do_reflection[5] = reflection_upper_z; + + do_stagger[0] = avoid_origin_x; + do_stagger[1] = avoid_origin_upper_x; + do_stagger[2] = avoid_origin_y; + do_stagger[3] = avoid_origin_upper_y; + do_stagger[4] = avoid_origin_z; + do_stagger[5] = avoid_origin_upper_z; + + + + /* Loop over all directions and faces */ + 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 we have the outer boundary of that face */ + if (cctkGH->cctk_bbox[2*dir+face]) + { + + /* Find parity */ + switch (ttype) + { + case SCALAR: + parity = +1; + break; + case VECTOR: + parity = dir == tcomponent ? -1 : +1; + break; + case TENSOR: + 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; + case 3: parity = +1; break; + case 4: parity = (dir == 1 || dir == 2) ? -1 : +1; break; + case 5: parity = +1; break; + default: assert (0); + } + break; + default: + assert (0); + } + + /* Find region extent */ + for (d=0; d<3; ++d) + { + lsh[d] = cctkGH->cctk_lsh[d]; + imin[d] = 0; + imax[d] = cctkGH->cctk_lsh[d]; + ioff[d] = 0; + } + if (face == 0) + { + imax[d] = cctkGH->cctk_nghostzones[d]; + ioff[d] = (+ cctkGH->cctk_nghostzones[d] + + (do_stagger[2*d+face] ? 0 : 1)); + } + else + { + imin[d] = cctkGH->cctk_lsh[d] - cctkGH->cctk_nghostzones[d]; + ioff[d] = (- cctkGH->cctk_nghostzones[d] + - (do_stagger[2*d+face] ? 0 : 1)); + } + + /* Copy region */ + switch (group.vartype) + { + +#ifdef CCTK_INT1 + case CCTK_VARIABLE_INT1: +#ifdef CCTK_INTER_PRECISION_1 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT1 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_INT2 + case CCTK_VARIABLE_INT2: +#ifdef CCTK_INTER_PRECISION_2 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT2 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_INT4 + case CCTK_VARIABLE_INT4: +#ifdef CCTK_INTER_PRECISION_4 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT4 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_INT8 + case CCTK_VARIABLE_INT8: +#ifdef CCTK_INTER_PRECISION_8 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT8 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL4 + case CCTK_VARIABLE_REAL4: +#ifdef CCTK_REAL_PRECISION_4 + case CCTK_VARIABLE_REAL: +#endif + copy_CCTK_REAL4 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL8 + case CCTK_VARIABLE_REAL8: +#ifdef CCTK_REAL_PRECISION_8 + case CCTK_VARIABLE_REAL: +#endif + copy_CCTK_REAL8 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL16 + case CCTK_VARIABLE_REAL16: +#ifdef CCTK_REAL_PRECISION_16 + case CCTK_VARIABLE_REAL: +#endif + copy_CCTK_REAL16 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL4 + case CCTK_VARIABLE_COMPLEX8: +#ifdef CCTK_COMPLEX_PRECISION_8 + case CCTK_VARIABLE_COMPLEX: +#endif + copy_CCTK_COMPLEX8 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL8 + case CCTK_VARIABLE_COMPLEX16: +#ifdef CCTK_COMPLEX_PRECISION_16 + case CCTK_VARIABLE_COMPLEX: +#endif + copy_CCTK_COMPLEX16 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL16 + case CCTK_VARIABLE_COMPLEX32: +#ifdef CCTK_COMPLEX_PRECISION_32 + case CCTK_VARIABLE_COMPLEX: +#endif + copy_CCTK_COMPLEX32 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + + default: + CCTK_WARN (0, "Unsupported variable type"); + } + + } /* if cctk_bbox */ + } /* if do_reflection */ + } /* for face */ + } /* for dir */ + + /* Success */ + return 0; +} + + + +void +ReflectionSymmetry_Apply (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + + int nvars; + CCTK_INT * restrict indices; + CCTK_INT * restrict faces; + CCTK_INT * restrict widths; + CCTK_INT * restrict tables; + int vi; + int dim; + int * restrict stencil; + int i; + int istat; + int ierr; + + if (! cctkGH) + { + CCTK_WARN (0, "Argument cctkGH is NULL"); + } + + nvars = Boundary_SelectedGVs (cctkGH, 0, NULL, NULL, NULL, NULL, NULL); + if (nvars < 0) + { + CCTK_WARN (0, "Internal error in Boundary_SelectedGVs"); + } + + indices = malloc (nvars * sizeof *indices); + if (! indices) + { + CCTK_WARN (0, "Out of memory"); + } + faces = malloc (nvars * sizeof *faces); + if (! faces) + { + CCTK_WARN (0, "Out of memory"); + } + widths = malloc (nvars * sizeof *widths); + if (! widths) + { + CCTK_WARN (0, "Out of memory"); + } + tables = malloc (nvars * sizeof *tables); + if (! tables) + { + CCTK_WARN (0, "Out of memory"); + } + + istat = Boundary_SelectedGVs + (cctkGH, nvars, indices, faces, widths, tables, 0); + if (istat != nvars) + { + CCTK_WARN (0, "Internal error in Boundary_SelectedGVs"); + } + + for (i=0; i= CCTK_NumVars()) + { + CCTK_WARN (0, "Illegal variable index"); + } + + if (widths[i] < 0) + { + CCTK_WARN (0, "Illegal boundary width"); + } + + dim = CCTK_GroupDimFromVarI (vi); + if (dim < 0) + { + CCTK_WARN (0, "Illegal dimension"); + } + + stencil = malloc (dim * sizeof *stencil); + if (! stencil) + { + CCTK_WARN (0, "Out of memory"); + } + ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi); + if (ierr) + { + CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI"); + } + + ierr = BndReflectVI (cctkGH, stencil, vi); + if (ierr) + { + CCTK_WARN (0, "Internal error in BndReflectVI"); + } + + free (stencil); + } + + free (indices); + free (faces); + free (widths); + free (tables); +} diff --git a/src/interpolate.c b/src/interpolate.c new file mode 100644 index 0000000..8591fd7 --- /dev/null +++ b/src/interpolate.c @@ -0,0 +1,359 @@ +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "reflection.h" + + + +static const char rcsid[] = "$Header$"; +CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_interpolate_c); + + + +CCTK_INT +ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, + CCTK_INT const N_dims, + CCTK_INT const local_interp_handle, + CCTK_INT const param_table_handle, + CCTK_INT const coord_system_handle, + CCTK_INT const N_interp_points, + CCTK_INT const interp_coords_type, + CCTK_POINTER_TO_CONST restrict const interp_coords[], + CCTK_INT const N_input_arrays, + CCTK_INT const input_array_indices[], + CCTK_INT const N_output_arrays, + CCTK_INT const output_array_types[], + CCTK_POINTER restrict const output_arrays[], + CCTK_INT const faces) +{ + cGH const * restrict const cctkGH = cctkGH_; + DECLARE_CCTK_PARAMETERS; + + int do_reflection[6]; + + CCTK_POINTER_TO_CONST restrict new_interp_coords[3]; + CCTK_INT newfaces; + + CCTK_INT * restrict operand_indices; + CCTK_INT * restrict operation_codes; + CCTK_INT * restrict output_array_indices; + + int dir; + + int iret; + + int m; + int n; + int d; + + int ierr; + + + + /* Get symmetry information */ + do_reflection[0] = reflection_x; + do_reflection[1] = reflection_upper_x; + do_reflection[2] = reflection_y; + do_reflection[3] = reflection_upper_y; + do_reflection[4] = reflection_z; + do_reflection[5] = reflection_upper_z; + + newfaces = faces; + for (d=0; d<6; ++d) + { + if (do_reflection[d]) + { + assert (newfaces & (1 << d)); + newfaces &= ~ (1 << d); + } + } + + + + /* Fold coordinates */ + assert (interp_coords_type == CCTK_VARIABLE_REAL); + for (dir=0; dir<3; ++dir) + { + assert (! do_reflection[2*dir+1]); + + if (do_reflection[2*dir]) + { + new_interp_coords[dir] + = malloc (N_interp_points * sizeof *new_interp_coords[dir]); + assert (new_interp_coords[d]); + + for (n=0; n=0 && operand_indices[m]=0 + && output_array_indices[m]=0 && vi=0 && gi0); + firstvar = CCTK_FirstVarIndexI(gi); + assert (firstvar>=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) { + 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, "u")) + { + /* vector */ + assert (numvars == 3); + ttype = VECTOR; + tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) { + /* symmetric tensor */ + assert (numvars == 6); + ttype = TENSOR; + 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; + default: + assert (0); + } + + + + /* 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; + 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; + } + } + + + + /* Are there negative parities? */ + needs_checking = 0; + for (dir=0; dir<3; ++dir) + { + check_dir[d] = do_reflection[2*dir] && parities[dir] < 0; + needs_checking |= check_dir[d]; + } + + + + /* Loop over all points and unfold */ + if (needs_checking) + { + for (n=0; n