interpolation.h File Reference

#include <petscpf.h>
#include <petscdmswarm.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <petsctime.h>
#include <petscdmcomposite.h>
#include <petscerror.h>
#include <petscsys.h>
#include "variables.h"
#include "ParticleSwarm.h"
#include "walkingsearch.h"
#include "grid.h"
#include "logging.h"
#include "io.h"
#include "setup.h"

Include dependency graph for interpolation.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	NUM_WEIGHTS   8
#define	InterpolateFieldFromCornerToCenter(field, centfield, user)
	Generic macro to call the appropriate interpolation function based on the field type.
#define	InterpolateFieldFromCenterToCorner(blockSize, centfield_ptr, corner_ptr, user_ctx)
	Macro to dispatch to the correct scalar or vector center-to-corner function based on a runtime block size variable.
#define	InterpolateCornerToFaceCenter(corner_arr, faceX_arr, faceY_arr, faceZ_arr, user_ctx)
	A type-generic macro that interpolates a field from corner nodes to all face centers.
#define	PieceWiseLinearInterpolation(fieldName, fieldPtr, i, j, k, outPtr)
	Macro that calls either the scalar or vector piecewise interpolation function based on the type of the fieldPtr parameter (3D array).
#define	TrilinearInterpolation(fieldName, fieldPtr, i, j, k, a1, a2, a3, outPtr)
	Macro that calls either the scalar or vector trilinear interpolation function based on the type of the fieldPtr parameter (3D array).

Functions
PetscErrorCode	TrilinearInterpolation_Scalar (const char *fieldName, PetscReal ***fieldScal, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, PetscReal *val)
	Computes the trilinear interpolated scalar at a given point.
PetscErrorCode	TrilinearInterpolation_Vector (const char *fieldName, Cmpnts ***fieldVec, PetscInt i, PetscInt j, PetscInt k, PetscReal a1, PetscReal a2, PetscReal a3, Cmpnts *vec)
	Computes the trilinear interpolated vector (e.g., velocity) at a given point.
PetscErrorCode	InterpolateEulerFieldToSwarm (UserCtx *user, Vec fieldLocal_cellCentered, const char *fieldName, const char *swarmOutFieldName)
	Interpolates a cell-centered field (scalar or vector) onto DMSwarm particles, using a robust, PETSc-idiomatic two-stage process.
PetscErrorCode	InterpolateAllFieldsToSwarm (UserCtx *user)
	Interpolates all relevant fields from the DMDA to the DMSwarm.
PetscErrorCode	InterpolateParticleVelocities (UserCtx *user)
	Interpolates particle velocities using trilinear interpolation.
PetscErrorCode	InterpolateFieldFromCornerToCenter_Scalar (PetscReal ***field, PetscReal ***centfield, UserCtx *user)
	Safely interpolate a scalar field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.
PetscErrorCode	InterpolateFieldFromCornerToCenter_Vector (Cmpnts ***field, Cmpnts ***centfield, UserCtx *user)
	Safely interpolate a vector field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.
PetscErrorCode	TestCornerToCenterInterpolation (UserCtx *user)
	Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Vector (Cmpnts ***centfield_arr, Cmpnts ***corner_arr, UserCtx *user)
	Interpolates a vector field from cell centers to corner nodes.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Scalar (PetscReal ***centfield_arr, PetscReal ***corner_arr, UserCtx *user)
	Interpolates a scalar field from cell centers to corner nodes.
PetscErrorCode	GetScatterTargetInfo (UserCtx *user, const char *particleFieldName, DM *targetDM, PetscInt *expected_dof)
	Determines the target Eulerian DM and expected DOF for scattering a given particle field.
PetscErrorCode	GetPersistentLocalVector (UserCtx *user, const char *fieldName, Vec *localVec)
	Retrieves the persistent local vector (e.g., lPsi, lUcat) for a given field name.
PetscErrorCode	AccumulateParticleField (DM swarm, const char *particleFieldName, DM gridSumDM, Vec gridSumVec)
	Accumulates a particle field (scalar or vector) into a target grid sum vector.
PetscErrorCode	NormalizeGridVectorByCount (DM countDM, Vec countVec, DM dataDM, Vec sumVec, Vec avgVec)
	Normalizes a grid vector of sums by a grid vector of counts to produce an average.
PetscErrorCode	ScatterParticleFieldToEulerField (UserCtx *user, const char *particleFieldName, Vec eulerFieldAverageVec)
	Scatters a particle field (scalar or vector) to the corresponding Eulerian field average.
PetscErrorCode	ScatterAllParticleFieldsToEulerFields (UserCtx *user)
	Scatters a predefined set of particle fields to their corresponding Eulerian fields.
PetscErrorCode	InterpolateCornerToFaceCenter_Scalar (PetscReal ***corner_arr, PetscReal ***faceX_arr, PetscReal ***faceY_arr, PetscReal ***faceZ_arr, UserCtx *user)
	Interpolates a scalar field from corner nodes to all face centers.
PetscErrorCode	InterpolateCornerToFaceCenter_Vector (Cmpnts ***corner_arr, Cmpnts ***faceX_arr, Cmpnts ***faceY_arr, Cmpnts ***faceZ_arr, UserCtx *user)
	Interpolates a vector field from corner nodes to all face centers.

Macro Definition Documentation

NUM_WEIGHTS

#define NUM_WEIGHTS   8

Definition at line 24 of file interpolation.h.

InterpolateFieldFromCornerToCenter

#define InterpolateFieldFromCornerToCenter	(		field,
			centfield,
			user
	)

Value:

  ( (void)sizeof(char[1 - 2*!!(!__builtin_types_compatible_p(typeof(field), typeof(centfield)))]),\
    _Generic((field),                                                                       \
      PetscReal ***: InterpolateFieldFromCornerToCenter_Scalar,                             \
      Cmpnts ***:    InterpolateFieldFromCornerToCenter_Vector                                \
    )(field, centfield, user) )

Generic macro to call the appropriate interpolation function based on the field type.

This macro will select either the scalar or vector interpolation function based on the type of the 'field' pointer. It also performs a compile-time check that 'centfield' is of the same type as 'field'. If the types are not the same, a compile-time error is produced.

Usage: InterpolateFieldFromCornerToCenter(field, centfield, user);

Definition at line 36 of file interpolation.h.

                   : InterpolateFieldFromCornerToCenter_Scalar,                             \
      Cmpnts ***:    InterpolateFieldFromCornerToCenter_Vector                                \
    )(field, centfield, user) )

InterpolateFieldFromCenterToCorner

#define InterpolateFieldFromCenterToCorner	(		blockSize,
			centfield_ptr,
			corner_ptr,
			user_ctx
	)

Value:

    ( (blockSize) == 1 ? \
        InterpolateFieldFromCenterToCorner_Scalar((PetscReal***)(centfield_ptr), (PetscReal***)(corner_ptr), (user_ctx)) : \
        InterpolateFieldFromCenterToCorner_Vector((Cmpnts***)(centfield_ptr), (Cmpnts***)(corner_ptr), (user_ctx)) \
    )

Macro to dispatch to the correct scalar or vector center-to-corner function based on a runtime block size variable.

This macro uses a ternary operator to inspect the runtime value of 'blockSize' and select the appropriate implementation. It expects the input and output pointers to be void* and handles casting them to the correct, strongly-typed pointers.

Parameters

blockSize	The runtime block size (1 for scalar, 3 for vector).
centfield_ptr	The input void* pointer to the 3D cell-centered data array.
corner_ptr	The output void* pointer to the 3D corner data array.
user_ctx	The UserCtx structure.

Definition at line 56 of file interpolation.h.

                                                                                                                         : \
        InterpolateFieldFromCenterToCorner_Vector((Cmpnts***)(centfield_ptr), (Cmpnts***)(corner_ptr), (user_ctx)) \
    )

InterpolateCornerToFaceCenter

#define InterpolateCornerToFaceCenter	(		corner_arr,
			faceX_arr,
			faceY_arr,
			faceZ_arr,
			user_ctx
	)

Value:

    _Generic((corner_arr),                                                                  \
        PetscReal***: InterpolateCornerToFaceCenter_Scalar,                                 \
        Cmpnts***:    InterpolateCornerToFaceCenter_Vector                                  \
    )(corner_arr, faceX_arr, faceY_arr, faceZ_arr, user_ctx)

A type-generic macro that interpolates a field from corner nodes to all face centers.

This macro uses C11's _Generic feature to dispatch to the appropriate underlying function based on the type of the input corner_arr.

• If corner_arr is of type PetscReal***, it calls InterpolateCornerToFaceCenter_Scalar.
• If corner_arr is of type Cmpnts***, it calls InterpolateCornerToFaceCenter_Vector.

Parameters

corner_arr	Ghosted node-centered array (global indexing). The type of this argument determines which function is called.
faceX_arr	Local array for X-faces.
faceY_arr	Local array for Y-faces.
faceZ_arr	Local array for Z-faces.
user_ctx	User context containing DMDA 'fda'.

Definition at line 79 of file interpolation.h.

                    : InterpolateCornerToFaceCenter_Scalar,                                 \
        Cmpnts***:    InterpolateCornerToFaceCenter_Vector                                  \
    )(corner_arr, faceX_arr, faceY_arr, faceZ_arr, user_ctx)

PieceWiseLinearInterpolation

#define PieceWiseLinearInterpolation	(		fieldName,
			fieldPtr,
			i,
			j,
			k,
			outPtr
	)

Value:

  _Generic((fieldPtr),                                                     \
    PetscReal ***: PieceWiseLinearInterpolation_Scalar,                    \
    Cmpnts      ***: PieceWiseLinearInterpolation_Vector                   \
  )(fieldName, fieldPtr, i, j, k, outPtr)

Macro that calls either the scalar or vector piecewise interpolation function based on the type of the fieldPtr parameter (3D array).

Usage example:

// For scalar: PetscReal ***fieldScal; PetscReal outVal; PieceWiseLinearInterpolation(fieldName, fieldScal, i, j, k, &outVal);

// For vector: Cmpnts ***fieldVec; Cmpnts vec; PieceWiseLinearInterpolation(fieldName, fieldVec, i, j, k, &vec);

Definition at line 101 of file interpolation.h.

                 : PieceWiseLinearInterpolation_Scalar,                    \
    Cmpnts      ***: PieceWiseLinearInterpolation_Vector                   \
  )(fieldName, fieldPtr, i, j, k, outPtr)

TrilinearInterpolation

#define TrilinearInterpolation	(		fieldName,
			fieldPtr,
			i,
			j,
			k,
			a1,
			a2,
			a3,
			outPtr
	)

Value:

  _Generic((fieldPtr),                                                                     \
    PetscReal ***: TrilinearInterpolation_Scalar,                                          \
    Cmpnts      ***: TrilinearInterpolation_Vector                                         \
  )(fieldName, fieldPtr, i, j, k, a1, a2, a3, outPtr)

Macro that calls either the scalar or vector trilinear interpolation function based on the type of the fieldPtr parameter (3D array).

Usage example:

PetscReal result; Cmpnts vec;

// For scalars: TrilinearInterpolation(fieldName, fieldScal, i, j, k, a1, a2, a3, &result);

// For vectors: TrilinearInterpolation(fieldName, fieldVec, i, j, k, a1, a2, a3, &vec);

Definition at line 122 of file interpolation.h.

                 : TrilinearInterpolation_Scalar,                                          \
    Cmpnts      ***: TrilinearInterpolation_Vector                                         \
  )(fieldName, fieldPtr, i, j, k, a1, a2, a3, outPtr)

Function Documentation

TrilinearInterpolation_Scalar()

PetscErrorCode TrilinearInterpolation_Scalar	(	const char *	fieldName,
		PetscReal ***	fieldScal,
		PetscInt	i,
		PetscInt	j,
		PetscInt	k,
		PetscReal	a1,
		PetscReal	a2,
		PetscReal	a3,
		PetscReal *	val
	)

Computes the trilinear interpolated scalar at a given point.

Parameters

[in]	fieldName	A string representing the name of the scalar field (e.g., "temperature").
[in]	fieldScal	3D array of the field from a DMDA (indexed as [k][j][i]), each cell a PetscReal.
[in]	i,j,k	Integral cell indices (the "lower" corner in each dimension).
[in]	a1,a2,a3	Normalized coordinates within the cell ([0,1] range).
[out]	val	Pointer to a PetscReal that will store the interpolated scalar.

This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same interface.

Computes the trilinear interpolated scalar at a given point.

Parameters

[in]	fieldName	A string representing the name of the scalar field (e.g., "temperature").
[in]	fieldScal	3D array of the field from a DMDA (indexed as [k][j][i]), each cell a PetscReal.
[in]	i,j,k	Integral cell indices (the "lower" corner in each dimension).
[in]	a1,a2,a3	Normalized coordinates within the cell ([0,1] range).
[out]	val	Pointer to a PetscReal that will store the Interpolateted scalar.

This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same PetscInterface.

Definition at line 661 of file interpolation.c.

{
    PetscFunctionBegin; // PETSc macro for error/stack tracing
 
    // Compute the 8 corner weights
    PetscReal wcorner[8];
    ComputeTrilinearWeights(a1, a2, a3, wcorner);
 
    // Offsets for cell corners
    PetscInt i1 = i + 1;
    PetscInt j1 = j + 1;
    PetscInt k1 = k + 1;
 
    // Initialize the output scalar
    PetscReal sum = 0.0;
 
    // Corner 0 => (i, j, k)
    sum += wcorner[0] * fieldScal[k ][j ][i ];
    // Corner 1 => (i+1, j, k)
    sum += wcorner[1] * fieldScal[k ][j ][i1];
    // Corner 2 => (i, j+1, k)
    sum += wcorner[2] * fieldScal[k ][j1][i ];
    // Corner 3 => (i+1, j+1, k)
    sum += wcorner[3] * fieldScal[k ][j1][i1];
    // Corner 4 => (i, j, k+1)
    sum += wcorner[4] * fieldScal[k1][j ][i ];
    // Corner 5 => (i+1, j, k+1)
    sum += wcorner[5] * fieldScal[k1][j ][i1];
    // Corner 6 => (i, j+1, k+1)
    sum += wcorner[6] * fieldScal[k1][j1][i ];
    // Corner 7 => (i+1, j+1, k+1)
    sum += wcorner[7] * fieldScal[k1][j1][i1];
 
    *val = sum;
 
    // Logging (optional)
    LOG_ALLOW(LOCAL, LOG_VERBOSE,
        "Field '%s' at (i=%d, j=%d, k=%d), "
        "a1=%.6f, a2=%.6f, a3=%.6f -> val=%.6f.\n",
        fieldName, i, j, k, a1, a2, a3, *val);
 
    // LOG_ALLOW_SYNC(GLOBAL, LOG_INFO,
    //    "TrilinearInterpolation_Scalar: Completed interpolation for field '%s' across local cells.\n",
      //    fieldName);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

TrilinearInterpolation_Vector()

PetscErrorCode TrilinearInterpolation_Vector	(	const char *	fieldName,
		Cmpnts ***	fieldVec,
		PetscInt	i,
		PetscInt	j,
		PetscInt	k,
		PetscReal	a1,
		PetscReal	a2,
		PetscReal	a3,
		Cmpnts *	vec
	)

Computes the trilinear interpolated vector (e.g., velocity) at a given point.

Parameters

[in]	fieldName	A string representing the name of the vector field (e.g., "velocity").
[in]	fieldVec	3D array of the field from a DMDA (indexed as [k][j][i]), each cell of type Cmpnts.
[in]	i,j,k	Integral cell indices (the "lower" corner in each dimension).
[in]	a1,a2,a3	Normalized coordinates within the cell ([0,1] range).
[out]	vec	Pointer to a Cmpnts struct that will store the interpolated vector (x, y, z).

This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same interface.

Computes the trilinear interpolated vector (e.g., velocity) at a given point.

Parameters

[in]	fieldName	A string representing the name of the vector field (e.g., "velocity").
[in]	fieldVec	3D array of the field from a DMDA (indexed as [k][j][i]), each cell of type Cmpnts.
[in]	i,j,k	Integral cell indices (the "lower" corner in each dimension).
[in]	a1,a2,a3	Normalized coordinates within the cell ([0,1] range).
[out]	vec	Pointer to a Cmpnts struct that will store the Interpolateted vector (x, y, z).

This function uses the standard 8-corner trilinear formula via ComputeTrilinearWeights(). If a different scheme is desired, implement a new function with the same PetscInterface.

Computes the trilinear Interpolateted vector at a given point, with partial weighting if corners go out of range.

If any of the 8 corners (i or i+1, j or j+1, k or k+1) is out of [0..mx), [0..my), [0..mz), that corner is skipped and the corresponding weight is omitted. The total is normalized by the sum of used weights, so we get a partial interpolation near boundaries.

Definition at line 741 of file interpolation.c.

{
  PetscFunctionBegin; // PETSc macro for error/stack tracing
 
  // Compute the 8 corner weights
  PetscErrorCode ierr;
  PetscReal wcorner[8];
  PetscMPIInt rank;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  
  LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] Computing Trilinear Weights.\n",rank);
  ComputeTrilinearWeights(a1, a2, a3, wcorner);
  
  LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d]  Trilinear weights computed for local cell %d,%d,%d.\n",rank,i,j,k);
 
  // For partial interpolation, we'll keep track of how many corners are valid
  // and how much sum of weights is used. Then we do a final normalization.
  PetscReal sumW = 0.0;
  Cmpnts    accum = {0.0, 0.0, 0.0};
 
  // The eight corner indices, with their weights:
  // corners: (i,j,k), (i+1,j,k), (i,j+1,k), (i+1,j+1,k), etc.
  // We store them in an array to iterate cleanly.
  const PetscInt cornerOffsets[8][3] = {
    {0, 0, 0},
    {1, 0, 0},
    {0, 1, 0},
    {1, 1, 0},
    {0, 0, 1},
    {1, 0, 1},
    {0, 1, 1},
    {1, 1, 1}
  };
 
  // Weighted partial sum
  for (PetscInt c = 0; c < 8; c++) {
    const PetscInt di = cornerOffsets[c][0];
    const PetscInt dj = cornerOffsets[c][1];
    const PetscInt dk = cornerOffsets[c][2];
    PetscInt iC = i + di;
    PetscInt jC = j + dj;
    PetscInt kC = k + dk;
 
    /*
    // skip if out of domain
    // (Assuming you know global domain is [0..mx), [0..my), [0..mz).)
    if (iC < 0 || iC >= (PetscInt)userGlobalMx ||
        jC < 0 || jC >= (PetscInt)userGlobalMy ||
        kC < 0 || kC >= (PetscInt)userGlobalMz)
    {
      // skip this corner
      continue;
    }
 
    */
 
    LOG_ALLOW(LOCAL,LOG_VERBOSE,"[Rank %d] %s[%d][%d][%d] = (%.4f,%.4f,%.4f).\n",rank,fieldName,kC,jC,iC,fieldVec[kC][jC][iC].x,fieldVec[kC][jC][iC].y,fieldVec[kC][jC][iC].z);
    
    // Otherwise, accumulate
    accum.x += wcorner[c] * fieldVec[kC][jC][iC].x;
    accum.y += wcorner[c] * fieldVec[kC][jC][iC].y;
    accum.z += wcorner[c] * fieldVec[kC][jC][iC].z;
    sumW += wcorner[c];
  }
 
  // If sumW=0 => out-of-range or zero weighting => set (0,0,0)
  if (sumW > 1.0e-14) {
    vec->x = accum.x / sumW;
    vec->y = accum.y / sumW;
    vec->z = accum.z / sumW;
  } else {
    vec->x = 0.0;  vec->y = 0.0;  vec->z = 0.0;
  }
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

InterpolateEulerFieldToSwarm()

PetscErrorCode InterpolateEulerFieldToSwarm	(	UserCtx *	user,
		Vec	fieldLocal_cellCentered,
		const char *	fieldName,
		const char *	swarmOutFieldName
	)

Interpolates a cell-centered field (scalar or vector) onto DMSwarm particles, using a robust, PETSc-idiomatic two-stage process.

This function first converts the cell-centered input data to corner-node data, storing this intermediate result in a PETSc Vec to correctly handle the communication of ghost-point information across parallel ranks. It then performs a final trilinear interpolation from the ghosted corner data to each particle's location.

Workflow:

1. Create temporary PETSc Vecs (cornerGlobal, cornerLocal) to manage the intermediate corner-node data, using the existing nodal DMDA (user->fda).
2. Call a dispatch macro that uses the runtime block size (bs) to select the correct underlying center-to-corner function, writing results into cornerGlobal.
3. Perform a ghost-point exchange (DMGlobalToLocal) to transfer the boundary data from cornerGlobal into the ghost regions of cornerLocal.
4. Loop over all local particles. For each particle: a. Convert its global cell index to a local index relative to the ghosted array. b. Check if the particle's interpolation stencil is fully contained within the owned+ghost region. If not, log a warning and set the result to zero. c. Perform the final trilinear interpolation using the ghosted cornerLocal data.
5. Restore all PETSc objects to prevent memory leaks.

Parameters

[in]	user	User context with DMDA, DMSwarm, etc.
[in]	fieldGlobal_cellCentered	Vec (from fda) with cell-centered data (e.g., Ucat).
[in]	fieldName	Human-readable field name for logging (e.g., "Ucat").
[in]	swarmOutFieldName	Name of the DMSwarm field where interpolation results go.

Returns

PetscErrorCode 0 on success.

This function first converts the cell-centered input data to corner-node data, storing this intermediate result in a PETSc Vec to correctly handle the communication of ghost-point information across parallel ranks. It then performs a final trilinear interpolation from the ghosted corner data to each particle's location.

Workflow:

1. Create temporary PETSc Vecs (cornerGlobal, cornerLocal) to manage the intermediate corner-node data, using the existing nodal DMDA (user->fda).
2. Call a dispatch macro that uses the runtime block size (bs) to select the correct underlying center-to-corner function, writing results into cornerGlobal.
3. Perform a ghost-point exchange (DMGlobalToLocal) to transfer the boundary data from cornerGlobal into the ghost regions of cornerLocal.
4. Loop over all local particles. For each particle: a. Convert its global cell index to a local index relative to the ghosted array. b. Check if the particle's interpolation stencil is fully contained within the owned+ghost region. If not, log a warning and set the result to zero. c. Perform the final trilinear interpolation using the ghosted cornerLocal data.
5. Restore all PETSc objects to prevent memory leaks.

Parameters

[in]	user	User context with DMDA, DMSwarm, etc.
[in]	fieldlocal_cellCentered	Vec (from fda) with cell-centered data (e.g., Ucat).
[in]	fieldName	Human-readable field name for logging (e.g., "Ucat").
[in]	swarmOutFieldName	Name of the DMSwarm field where interpolation results go.

Returns

PetscErrorCode 0 on success.

Definition at line 997 of file interpolation.c.

{
    PetscErrorCode ierr;
    DM             fda = user->fda;
    DM             swarm = user->swarm;
    PetscInt       bs;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
 
    // Generic pointers to the raw data arrays
    void *cellCenterPtr_read;
    void *cornerPtr_read_with_ghosts;
 
    // Swarm-related pointers
    PetscInt  *cellIDs = NULL;
    PetscReal *weights = NULL;
    PetscInt64  *pids = NULL; // For logging particle IDs in warnings
    void      *swarmOut = NULL;
    PetscReal *pos = NULL;
    PetscInt  *status = NULL;
    PetscInt   nLocal;
 
    PetscFunctionBegin;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(fda, &info); CHKERRQ(ierr);
    ierr = VecGetBlockSize(fieldLocal_cellCentered, &bs); CHKERRQ(ierr);
    if (bs != 1 && bs != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "BlockSize must be 1 or 3.");
 
    LOG_ALLOW(GLOBAL, LOG_INFO, "Starting for field '%s'... \n", fieldName);
    // Select the appropriate DM for intermediate corner data based on block size
    DM dm_corner = (bs == 3) ? user->fda : user->da;
 
    /* STAGE 1: Center-to-Corner Calculation with Communication */
 
    // (A) Create or reuse the corner vectors.
    if(user->CellFieldAtCorner){
        PetscInt cached_bs;
        ierr = VecGetBlockSize(user->CellFieldAtCorner, &cached_bs); CHKERRQ(ierr);
        if(cached_bs != bs){
                LOG_ALLOW(LOCAL, LOG_WARNING, "[Rank %d] Block size mismatch for existing corner vector and current field %s(Cached: %d, New: %d). Replacing.\n", rank, fieldName, cached_bs, bs);
            ierr = VecDestroy(&user->CellFieldAtCorner); CHKERRQ(ierr);
            ierr = VecDestroy(&user->lCellFieldAtCorner); CHKERRQ(ierr);
            user->CellFieldAtCorner = NULL;
            user->lCellFieldAtCorner = NULL;    
        }
    }
 
    if(user->CellFieldAtCorner == NULL){
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new global corner vector for field '%s'.\n", rank, fieldName);
        ierr = DMCreateGlobalVector(dm_corner, &user->CellFieldAtCorner); CHKERRQ(ierr);
    }else{
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing global corner vector for field '%s'.\n", rank, fieldName);
    }
 
    ierr = VecSet(user->CellFieldAtCorner,0.0); CHKERRQ(ierr);
    Vec cornerGlobal = user->CellFieldAtCorner;
 
    if(user->lCellFieldAtCorner == NULL){
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new local corner vector for field '%s'.\n", rank, fieldName);
        ierr = DMCreateLocalVector(dm_corner, &user->lCellFieldAtCorner); CHKERRQ(ierr);
    }else{
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Reusing existing local corner vector for field '%s'.\n", rank, fieldName);
    }
 
    ierr = VecSet(user->lCellFieldAtCorner,0.0); CHKERRQ(ierr);
    Vec cornerLocal = user->lCellFieldAtCorner;
 
    // (B) Get a read-only array from the input cell-centered vector
    ierr = DMDAVecGetArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
 
    PetscInt size; 
    ierr = VecGetSize(cornerGlobal, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
    
    size = 0;
 
    ierr = VecGetSize(cornerLocal, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
 
    size = 0;
    
    ierr = VecGetSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
    
    PetscInt xs,ys,zs,gxs,gys,gzs;
    
    ierr  = DMDAGetCorners(dm_corner,&xs,&ys,&zs,NULL,NULL,NULL); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] DMDAGetCorners for field '%s': xs=%d, ys=%d, zs=%d.\n", rank, fieldName, (PetscInt)xs, (PetscInt)ys, (PetscInt)zs);
 
    ierr  = DMDAGetGhostCorners(dm_corner,&gxs,&gys,&gzs,NULL,NULL,NULL); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] DMDAGetGhostCorners for field '%s': gxs=%d, gys=%d, gzs=%d.\n", rank, fieldName, (PetscInt)gxs, (PetscInt)gys, (PetscInt)gzs);
    
    // DEBUG: Inspect Ucat ghost cells before interpolation
    /*
    if (bs == 3) {
        Cmpnts ***ucat_array = (Cmpnts***)cellCenterPtr_read;
        DMDALocalInfo info_debug;
        ierr = DMDAGetLocalInfo(fda, &info_debug); CHKERRQ(ierr);
 
        // Only print on rank 0 to avoid clutter
        if (rank == 0) {
            PetscPrintf(PETSC_COMM_SELF, "\nDEBUG: Inspecting Ucat Ghost Cells...\n");
            PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n");
 
            // --- Check the MINIMUM-SIDE (-Xi) ---
            // The first physical cell is at index info.xs (local index).
            // The first ghost cell is at info.xs - 1.
            PetscInt i_first_phys = info_debug.xs;
            PetscInt i_first_ghost = info_debug.xs - 1;
            PetscInt j_mid = info_debug.ys + info_debug.ym / 2;
            PetscInt k_mid = info_debug.zs + info_debug.zm / 2;
 
            PetscPrintf(PETSC_COMM_SELF, "MIN-SIDE (-Xi):\n");
            PetscPrintf(PETSC_COMM_SELF, "  Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_first_ghost,
                        ucat_array[k_mid][j_mid][i_first_ghost].x,
                        ucat_array[k_mid][j_mid][i_first_ghost].y,
                        ucat_array[k_mid][j_mid][i_first_ghost].z);
            PetscPrintf(PETSC_COMM_SELF, "  First Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_first_phys,
                        ucat_array[k_mid][j_mid][i_first_phys].x,
                        ucat_array[k_mid][j_mid][i_first_phys].y,
                        ucat_array[k_mid][j_mid][i_first_phys].z);
 
 
            // --- Check the MAXIMUM-SIDE (+Xi) ---
            // The last physical cell is at index info.xs + info.xm - 1.
            // The first ghost cell on the max side is at info.xs + info.xm.
            PetscInt i_last_phys = info_debug.xs + info_debug.xm - 1;
            PetscInt i_last_ghost = info_debug.xs + info_debug.xm;
 
            PetscPrintf(PETSC_COMM_SELF, "MAX-SIDE (+Xi):\n");
            PetscPrintf(PETSC_COMM_SELF, "  Last Phys Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_last_phys,
                        ucat_array[k_mid][j_mid][i_last_phys].x,
                        ucat_array[k_mid][j_mid][i_last_phys].y,
                        ucat_array[k_mid][j_mid][i_last_phys].z);
            PetscPrintf(PETSC_COMM_SELF, "  Ghost Cell Ucat[%d][%d][%d] = (% .6f, % .6f, % .6f)\n",
                        k_mid, j_mid, i_last_ghost,
                        ucat_array[k_mid][j_mid][i_last_ghost].x,
                        ucat_array[k_mid][j_mid][i_last_ghost].y,
                        ucat_array[k_mid][j_mid][i_last_ghost].z);
            PetscPrintf(PETSC_COMM_SELF, "--------------------------------------------------\n\n");
        }
    }
    */
    // DEBUG   
 
    // (C) Perform the center-to-corner interpolation directly into the global corner vector
    void *cornerPtr_write = NULL;
    ierr = DMDAVecGetArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Starting center-to-corner interpolation for '%s'.\n", rank, fieldName);
        
    // SINGLE, CLEAN CALL SITE: The macro handles the runtime dispatch based on 'bs'.
    ierr = InterpolateFieldFromCenterToCorner(bs, cellCenterPtr_read, cornerPtr_write, user); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Finished center-to-corner interpolation for '%s'.\n", rank, fieldName);
    ierr = DMDAVecRestoreArray(dm_corner, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
    
    ierr = DMDAVecRestoreArrayRead(dm_corner, fieldLocal_cellCentered, &cellCenterPtr_read); CHKERRQ(ierr);
 
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    //////////// DEBUG
    /*
        // Log a synchronized header message from all ranks.
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerGlobal BEFORE ghost exchange...\n");
 
    // VecView prints to a PETSc viewer. PETSC_VIEWER_STDOUT_WORLD is a global, synchronized viewer.
    // We only need to call it if logging is active for this function.
    if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
        ierr = VecView(cornerGlobal, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
    
    // Log a synchronized footer message.
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerGlobal.\n");
    */
    ////////// DEBUG
 
    // ierr = PetscBarrier((PetscObject)cornerGlobal); CHKERRQ(ierr);
    
    // (D) CRITICAL STEP: Communicate the newly computed corner data to fill ghost regions
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Beginning ghost exchange for corner data...\n", rank);
 
    ierr = DMGlobalToLocalBegin(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(dm_corner, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
    
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Ghost exchange for corner data complete.\n", rank);
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_VERBOSE, "getting array from %s on all ranks.\n", dm_corner == user->fda ? "fda" : "da");
    if (is_function_allowed(__func__) && (int)(LOG_VERBOSE) <= (int)get_log_level()) {
        // DEBUG: Inspect cornerLocal after ghost exchange
        
        ierr = LOG_FIELD_ANATOMY(user,"CornerField", "After Corner Velocity Interpolated"); CHKERRQ(ierr);
 
        // DEBUG
        //ierr = DMView(user->fda, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    }
    
    /////// DEBUG
    /*
        // Log a synchronized header message from all ranks.
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Dumping cornerLocal AFTER ghost exchange...\n");
 
    // Here, we want each rank to print its own local vector.
    // PETSC_VIEWER_STDOUT_SELF is the correct tool for this. The LOG_ALLOW_SYNC
    // wrapper will ensure the output from each rank is printed sequentially.
    if (is_function_allowed(__func__) && (int)(LOG_DEBUG) <= (int)get_log_level()) {
        PetscMPIInt rank_d;
        ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank_d); CHKERRQ(ierr);
        
        // Use a synchronized print to create a clean header for each rank's output
        PetscSynchronizedPrintf(PETSC_COMM_WORLD, "\n--- cornerLocal (Rank %d) ---\n", rank_d);
        PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
 
        // Print the local vector's contents for this rank
        ierr = VecView(cornerLocal, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
 
        // Use a synchronized print to create a clean footer
        PetscSynchronizedPrintf(PETSC_COMM_WORLD, "--- End cornerLocal (Rank %d) ---\n", rank_d);
        PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG, "DEBUG: Finished dumping cornerLocal.\n");
    */
    ///// DEBUG
    
    // STAGE 2: Particle Interpolation using Ghosted Corner Data */
 
    // (E) Get the local, GHOSTED array of corner data for the final interpolation step
    ierr = DMDAVecGetArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
 
    // (F) Retrieve swarm fields for the particle loop
    ierr = DMSwarmGetLocalSize(swarm, &nLocal); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
    ierr = DMSwarmGetField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
 
    LOG_ALLOW(LOCAL,LOG_TRACE," Rank %d holds data upto & including %d,%d,%d.\n",rank,info.gxs + info.gxm,info.gys+info.gym,info.gzs+info.gzm);
 
    // (G) Loop over each local particle
    for (PetscInt p = 0; p < nLocal; p++) {
        
        Particle particle;
 
        ierr = UnpackSwarmFields(p,pids,weights,pos,cellIDs,NULL,status,NULL,NULL,NULL,&particle); CHKERRQ(ierr);
        
        LOG_ALLOW(LOCAL, LOG_VERBOSE,
          "[Rank %d] Particle PID %lld: global cell=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
          rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2],
          particle.weights.x, particle.weights.y, particle.weights.z);
        // Safety check: Ensure the entire 8-node stencil is within the valid memory
        // region (owned + ghosts) of the local array.
 
        if (particle.cell[0] < info.gxs || particle.cell[0] >= info.gxs + info.gxm - 1 ||
            particle.cell[1] < info.gys || particle.cell[1] >= info.gys + info.gym - 1 ||
            particle.cell[2] < info.gzs || particle.cell[2] >= info.gzs + info.gzm - 1)
        {
            LOG_ALLOW(LOCAL, LOG_WARNING,
                      "[Rank %d] Particle PID %lld in global cell (%d,%d,%d) is in an un-interpolatable region (requires ghosts of ghosts or is out of bounds). Zeroing field '%s'.\n",
                      rank, (long long)particle.PID, particle.cell[0], particle.cell[1], particle.cell[2], fieldName);
            if (bs == 3) {
                ((PetscReal*)swarmOut)[3*p + 0] = 0.0;
                ((PetscReal*)swarmOut)[3*p + 1] = 0.0;
                ((PetscReal*)swarmOut)[3*p + 2] = 0.0;
            } else {
                 ((PetscReal*)swarmOut)[p] = 0.0;
            }
            continue;
        }
 
        ierr = InterpolateEulerFieldToSwarmForParticle(
                  fieldName,
                  cornerPtr_read_with_ghosts,
                  &particle,
                  swarmOut, p, bs);
        CHKERRQ(ierr);
    }
 
    // (H) Restore all retrieved arrays and vectors */
    ierr = DMDAVecRestoreArrayRead(dm_corner, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
    
    // (I) Restore swarm fields
    ierr = DMSwarmRestoreField(swarm, swarmOutFieldName, NULL, NULL, &swarmOut); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "weight", NULL, NULL, (void**)&weights); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_pid", NULL, NULL, (void**)&pids); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_CellID", NULL, NULL, (void**)&cellIDs); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "position", NULL, NULL, (void**)&pos); CHKERRQ(ierr);
    ierr = DMSwarmRestoreField(swarm, "DMSwarm_location_status", NULL, NULL, (void**)&status); CHKERRQ(ierr);
    LOG_ALLOW(GLOBAL, LOG_INFO, "Finished for field '%s'.\n", fieldName);
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InterpolateAllFieldsToSwarm()

PetscErrorCode InterpolateAllFieldsToSwarm

(

UserCtx *

user

)

Interpolates all relevant fields from the DMDA to the DMSwarm.

Currently, it interpolates:

• user->Ucat (vector field) into the DMSwarm field "swarmVelocity".

To add more fields, duplicate the call to InterpolateOneFieldOverSwarm and provide:

• The global Vec for that field (e.g. user->Tcat for temperature),
• A human-readable field name (for logging),
• A DMSwarm output field name (e.g. "swarmTemperature").

Parameters

[in,out]

user

Pointer to a UserCtx containing: user->da (DM for the grid), user->swarm (DMSwarm for particles), user->Ucat (Vec for the vector field), possibly more fields like user->Tcat, user->Pcat, etc.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Currently, it Interpolatetes:

• user->Ucat (vector field) into the DMSwarm field "swarmVelocity".

To add more fields, duplicate the call to InterpolateOneFieldOverSwarm and provide:

• The global Vec for that field (e.g. user->Tcat for temperature),
• A human-readable field name (for logging),
• A DMSwarm output field name (e.g. "swarmTemperature").

Parameters

[in,out]

user

Pointer to a UserCtx containing: user->da (DM for the grid), user->swarm (DMSwarm for particles),. user->Ucat (Vec for the vector field), possibly more fields like user->Tcat, user->Pcat, etc.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

Definition at line 1321 of file interpolation.c.

{
  PetscErrorCode ierr;
  PetscMPIInt rank;
  PetscFunctionBegin;
 
  PROFILE_FUNCTION_BEGIN;
 
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank); CHKERRQ(ierr);
 
  /* 
     1) Interpolate the 'velocity' field (user->Ucat) into DMSwarm's "swarmVelocity".
        - The function InterpolateOneFieldOverSwarm() loops over all local particles,
          retrieves iCell, jCell, kCell, (a1,a2,a3) from the DMSwarm,
          and calls the appropriate trilinear interpolation routine.
 
        - "velocity" => just a label used in logging or debugging.
        - "swarmVelocity" => the DMSwarm field name where we store the result
          (assume dof=3 if user->Ucat has blockSize=3).
  */
 
  LOG_ALLOW(LOCAL, LOG_DEBUG,
         " Interpolation of ucat to velocity begins on rank %d.\n",rank);
  // Make sure to pass the *LOCAL* Vector to the function below! 
  ierr = InterpolateEulerFieldToSwarm(user, user->lUcat, 
                                      "Ucat", 
                                      "velocity"); CHKERRQ(ierr);
  ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivity,
                                      "Diffusivity", "Diffusivity"); CHKERRQ(ierr);
  
  ierr = InterpolateEulerFieldToSwarm(user,user->lDiffusivityGradient,
                                      "DiffusivityGradient", "DiffusivityGradient"); CHKERRQ(ierr);
  /* Add fields as necessary here*/
                                      
  ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
  LOG_ALLOW(LOCAL, LOG_INFO,
        "[rank %d]Completed Interpolateting all fields to the swarm.\n",rank);
 
  PROFILE_FUNCTION_END;
 
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

InterpolateParticleVelocities()

PetscErrorCode InterpolateParticleVelocities

(

UserCtx *

user

)

Interpolates particle velocities using trilinear interpolation.

Parameters

[in]

user

Pointer to the user-defined context containing grid and swarm information.

Returns

PetscErrorCode Returns 0 on success, non-zero on failure.

InterpolateFieldFromCornerToCenter_Scalar()

PetscErrorCode InterpolateFieldFromCornerToCenter_Scalar	(	PetscReal ***	field_arr,
		PetscReal ***	centfield_arr,
		UserCtx *	user
	)

Safely interpolate a scalar field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.

For each cell center in the physical region of the cell-centered DM (fda), this function averages the up to 8 surrounding scalar values from the coordinate DM (da). On boundaries, where fewer corners are available, a partial average is computed.

The coordinate DM (da) is built on corners (IM+1 x JM+1 x KM+1) while the cell-centered DM (fda) covers the physical cells (IM x JM x KM). Index offsets are adjusted via DMDAGetLocalInfo.

Parameters

[in]	field	3D array of corner-based scalar data (from user->da).
[out]	centfield	3D array for the interpolated cell-center scalar data (for user->fda).
[in]	user	User context containing: da : DM for the coordinate (corner) data. fda : DM for the cell-centered data.

Returns

PetscErrorCode 0 on success.

Safely interpolate a scalar field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.

For each interior computational cell (k,j,i), this function calculates the cell-centered value as the arithmetic average of the values at its 8 surrounding corner nodes.

This function adheres to the convention that the valid computational domain for cell-centered quantities excludes the first and last indices in each direction (e.g., cells at i=0 and i=IM). The loop bounds are constructed to match this.

Parameters

[in]	field_arr	Input: Ghosted 3D array of corner-node data, accessed via GLOBAL indices.
[out]	centfield_arr	Output: 3D array for interpolated cell-center values, accessed via GLOBAL indices.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Definition at line 111 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAGetLocalInfo(user->da, &info); CHKERRQ(ierr);
    
    // Get local and global grid dimensions
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    // Determine loop bounds to compute for interior cells only, matching the code's convention.
    // Start at index 1 if the process owns the global boundary at index 0.
    PetscInt is = (xs == 0) ? 1 : xs;
    PetscInt js = (ys == 0) ? 1 : ys;
    PetscInt ks = (zs == 0) ? 1 : zs;
 
    // Stop one cell short if the process owns the global boundary at the max index.
    PetscInt ie = (xe == mx) ? xe - 1 : xe;
    PetscInt je = (ye == my) ? ye - 1 : ye;
    PetscInt ke = (ze == mz) ? ze - 1 : ze;
 
    // Loop over the locally owned INTERIOR cells.
    for (PetscInt k = ks; k < ke; k++) {
        for (PetscInt j = js; j < je; j++) {
            for (PetscInt i = is; i < ie; i++) {
                // Calculate cell center value as the average of its 8 corner nodes
                centfield_arr[k][j][i] = 0.125 * (field_arr[k][j][i] + field_arr[k][j-1][i] +
                                                  field_arr[k-1][j][i] + field_arr[k-1][j-1][i] +
                                                  field_arr[k][j][i-1] + field_arr[k][j-1][i-1] +
                                                  field_arr[k-1][j][i-1] + field_arr[k-1][j-1][i-1]);
            }
        }
    }
 
    PetscFunctionReturn(0);
}

InterpolateFieldFromCornerToCenter_Vector()

PetscErrorCode InterpolateFieldFromCornerToCenter_Vector	(	Cmpnts ***	field_arr,
		Cmpnts ***	centfield_arr,
		UserCtx *	user
	)

Safely interpolate a vector field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.

For each cell center in the physical region of the cell-centered DM (fda), this function averages the 8 surrounding corner values from the coordinate DM (da). The coordinate DM (da) is built on corners (IM+1 x JM+1 x KM+1) while the cell-centered DM (fda) covers the physical cells (IM x JM x KM). Index offsets are adjusted using DMDAGetLocalInfo.

Parameters

[in]	field	3D array of corner-based vector data (from user->da).
[out]	centfield	3D array for the interpolated cell-center vector data (for user->fda).
[in]	user	User context containing: da : DM for the coordinate (corner) data. fda : DM for the cell-centered data.

Returns

PetscErrorCode 0 on success.

Safely interpolate a vector field from corner nodes (from the coordinate DM) to cell centers (from the cell-centered DM) using the provided UserCtx.

For each interior computational cell (k,j,i), this function calculates the cell-centered value as the arithmetic average of the values at its 8 surrounding corner nodes.

This function adheres to the convention that the valid computational domain for cell-centered quantities excludes the first and last indices in each direction (e.g., cells at i=0 and i=IM). The loop bounds are constructed to match this, making it a direct counterpart to functions like ComputeCellCentersAndSpacing.

Parameters

[in]	field_arr	Input: Ghosted 3D array of corner-node data, accessed via GLOBAL indices.
[out]	centfield_arr	Output: 3D array for interpolated cell-center values, accessed via GLOBAL indices.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Definition at line 38 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
 
    PetscFunctionBeginUser;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Get local and global grid dimensions
    PetscInt xs = info.xs, xe = info.xs + info.xm;
    PetscInt ys = info.ys, ye = info.ys + info.ym;
    PetscInt zs = info.zs, ze = info.zs + info.zm;
    PetscInt mx = info.mx, my = info.my, mz = info.mz;
 
    // Determine loop bounds to compute for interior cells only, matching the code's convention.
    // Start at index 1 if the process owns the global boundary at index 0.
    PetscInt is = (xs == 0) ? 1 : xs;
    PetscInt js = (ys == 0) ? 1 : ys;
    PetscInt ks = (zs == 0) ? 1 : zs;
 
    // Stop one cell short if the process owns the global boundary at the max index.
    PetscInt ie = (xe == mx) ? xe - 1 : xe;
    PetscInt je = (ye == my) ? ye - 1 : ye;
    PetscInt ke = (ze == mz) ? ze - 1 : ze;
 
    // Loop over the locally owned INTERIOR cells.
    for (PetscInt k = ks; k < ke; k++) {
        for (PetscInt j = js; j < je; j++) {
            for (PetscInt i = is; i < ie; i++) {
                // Calculate cell center value as the average of its 8 corner nodes
                centfield_arr[k][j][i].x = 0.125 * (field_arr[k][j][i].x + field_arr[k][j-1][i].x +
                                                    field_arr[k-1][j][i].x + field_arr[k-1][j-1][i].x +
                                                    field_arr[k][j][i-1].x + field_arr[k][j-1][i-1].x +
                                                    field_arr[k-1][j][i-1].x + field_arr[k-1][j-1][i-1].x);
 
                centfield_arr[k][j][i].y = 0.125 * (field_arr[k][j][i].y + field_arr[k][j-1][i].y +
                                                    field_arr[k-1][j][i].y + field_arr[k-1][j-1][i].y +
                                                    field_arr[k][j][i-1].y + field_arr[k][j-1][i-1].y +
                                                    field_arr[k-1][j][i-1].y + field_arr[k-1][j-1][i-1].y);
 
                centfield_arr[k][j][i].z = 0.125 * (field_arr[k][j][i].z + field_arr[k][j-1][i].z +
                                                    field_arr[k-1][j][i].z + field_arr[k-1][j-1][i].z +
                                                    field_arr[k][j][i-1].z + field_arr[k][j-1][i-1].z +
                                                    field_arr[k-1][j][i-1].z + field_arr[k-1][j-1][i-1].z);
            }
        }
    }
 
    PetscFunctionReturn(0);
}

TestCornerToCenterInterpolation()

PetscErrorCode TestCornerToCenterInterpolation

(

UserCtx *

user

)

Tests the InterpolateFieldFromCornerToCenter function by reproducing the Cent vector.

This function serves as a unit test. It performs the following steps:

1. Takes the corner-centered nodal coordinates (from DMGetCoordinatesLocal) as input.
2. Uses the InterpolateFieldFromCornerToCenter macro to interpolate these coordinates to the cell centers, storing the result in a new temporary vector.
3. Compares this new vector with the user->Cent vector, which is assumed to have been computed by ComputeCellCentersAndSpacing and serves as the ground truth.
4. A 2-norm of the difference is computed. If it is below a small tolerance, the test passes.

Note

This function should be called immediately after ComputeCellCentersAndSpacing has been successfully executed.

Parameters

user	The UserCtx for a specific grid level.

Returns

PetscErrorCode 0 on success, or a PETSc error code on failure.

Definition at line 175 of file interpolation.c.

{
    PetscErrorCode ierr;
    Vec            lCoords, TestCent;
    Cmpnts         ***coor_arr, ***test_cent_arr;
    PetscReal      diff_norm;
 
    PetscFunctionBeginUser;
 
    // 1. Create a temporary vector to hold the result of our interpolation.
    //    It must have the same layout and size as the ground-truth user->Cent vector.
    ierr = VecDuplicate(user->Cent, &TestCent); CHKERRQ(ierr);
 
    // 2. Get the input (corner coordinates) and output (our test vector) arrays.
    ierr = DMGetCoordinatesLocal(user->da, &lCoords); CHKERRQ(ierr);
    ierr = DMDAVecGetArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecGetArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
 
    // 3. Call the generic interpolation macro.
    //    The macro will see that `coor_arr` is of type `Cmpnts***` and correctly
    //    call the `InterpolateFieldFromCornerToCenter_Vector` function.
    ierr = InterpolateFieldFromCornerToCenter(coor_arr, test_cent_arr, user); CHKERRQ(ierr);
 
    // 4. Restore the arrays.
    ierr = DMDAVecRestoreArrayRead(user->fda, lCoords, &coor_arr); CHKERRQ(ierr);
    ierr = DMDAVecRestoreArray(user->fda, TestCent, &test_cent_arr); CHKERRQ(ierr);
 
    // 5. IMPORTANT: Assemble the vector so its values are communicated across processors
    //    and it's ready for global operations like VecNorm.
    ierr = VecAssemblyBegin(TestCent); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(TestCent); CHKERRQ(ierr);
 
    // 6. Compare the result with the ground truth.
    //    We compute TestCent = -1.0 * user->Cent + 1.0 * TestCent.
    //    This calculates the difference vector: TestCent - user->Cent.
    ierr = VecAXPY(TestCent, -1.0, user->Cent); CHKERRQ(ierr);
 
    //    Now, compute the L2 norm of the difference vector. If the functions are
    //    identical, the norm should be zero (or very close due to floating point).
    ierr = VecNorm(TestCent, NORM_2, &diff_norm); CHKERRQ(ierr);
 
    // 7. Report the result and clean up.
    if (diff_norm < 1.0e-12) {
        LOG_ALLOW(GLOBAL,LOG_DEBUG,"[SUCCESS] Test passed. Norm of difference is %g.\n", (double)diff_norm);
    } else {
        LOG_ALLOW(GLOBAL,LOG_DEBUG, "[FAILURE] Test failed. Norm of difference is %g.\n", (double)diff_norm);
    }
 
    ierr = VecDestroy(&TestCent); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

InterpolateFieldFromCenterToCorner_Vector()

PetscErrorCode InterpolateFieldFromCenterToCorner_Vector	(	Cmpnts ***	centfield_arr,
		Cmpnts ***	corner_arr,
		UserCtx *	user
	)

Interpolates a vector field from cell centers to corner nodes.

This version is adapted to write directly into a ghosted local array obtained from DMDAVecGetArray(), which allows using GLOBAL indices for writing to the OWNED portion of the array.

Parameters

[in]	centfield_arr	Input: 3D array (ghosted) of cell-centered data, accessed via GLOBAL indices.
[out]	corner_arr	Output: 3D array (ghosted) where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Parameters

[in]	centfield_arr	Input: 3D array (ghosted) of cell-centered data, accessed via GLOBAL indices.
[out]	corner_arr	Output: 3D array where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Definition at line 241 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
    PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
    PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
 
    PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
    PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
    PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
 
 
    // Global grid dimensions (used for valid cell check)
    PetscInt IM = info.mx - 1; // Total nodes in i-direction
    PetscInt JM = info.my - 1; // Total nodes in j-direction
    PetscInt KM = info.mz - 1; // Total nodes in k-direction
 
    // Ghosted array starting indices (from the same info struct)
    PetscInt gxs = info.gxs;
    PetscInt gys = info.gys;
    PetscInt gzs = info.gzs;
 
    PetscInt xs = info.xs;
    PetscInt ys = info.ys;
    PetscInt zs = info.zs;
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
        "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
        rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
 
    // Loop over the GLOBAL indices of the NODES owned by this processor
    for (PetscInt k = zs_node; k < ze_node; k++) {
        for (PetscInt j = ys_node; j < ye_node; j++) {
            for (PetscInt i = xs_node; i < xe_node; i++) {
                Cmpnts sum = {0.0, 0.0, 0.0};
                PetscInt count = 0;
 
                // DEBUG 1 TEST
                /*
                if(rank == 1 && i == 24 && j == 12 && k == 49){
                    PetscInt i_cell_A = i - 1; 
                    PetscInt i_cell_B = i;
 
                    Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
                    Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
 
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Read { Valid Cell }Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f) and {Out-of-Bounds Cell}Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
                             k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
 
                }
                */
                
                // Skip processing the unused last node in each dimension.
                if(i >= IM || j >= JM || k >= KM){
                    continue;
                }
                // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
                // The index offsets correspond to the relative position of the indices(shifted) that represent cell-centered field values of cells that share the node N(k,j,i) as a corner
                for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
                    for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
                        for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
 
                            // These are still GLOBAL cell indices
                            PetscInt global_cell_k = k + dk_offset;
                            PetscInt global_cell_j = j + dj_offset;
                            PetscInt global_cell_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index
                            if (global_cell_i >= 0 && global_cell_i < nCellsX &&
                                global_cell_j >= 0 && global_cell_j < nCellsY &&
                                global_cell_k >= 0 && global_cell_k < nCellsZ)
                            {
                                    Cmpnts cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
 
                                    LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f, %.2f, %.2f)\n",
                                        rank,global_cell_k,global_cell_j,global_cell_i,cell_val.x, cell_val.y, cell_val.z);
                    
                                    sum.x += cell_val.x;
                                    sum.y += cell_val.y;
                                    sum.z += cell_val.z;
                                    count++;                                  
                            }
                        }
                    }
                }
 
                PetscInt i_global_write = i; // Global index in GLOBAL array.
                PetscInt j_global_write = j;
                PetscInt k_global_write = k;
                
                // We write directly into the array using the global loop indices.
                if (count > 0) {
                    corner_arr[k_global_write][j_global_write][i_global_write].x = sum.x / (PetscReal)count;
                    corner_arr[k_global_write][j_global_write][i_global_write].y = sum.y / (PetscReal)count;
                    corner_arr[k_global_write][j_global_write][i_global_write].z = sum.z / (PetscReal)count;
                } else {
                    // This case should ideally not happen for a valid owned node, but as a failsafe:
                    corner_arr[k_global_write][j_global_write][i_global_write] = (Cmpnts){0.0, 0.0, 0.0};
                }
               
                // DEBUG 2
                /*
                if(rank == 1){
                    if(i == 11 && j == 11 && k == 49){
                        Cmpnts ucat_node = corner_arr[k][j][i];
                        PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                                 rank, k, j, i,
                                 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                    }
                }
 
                if(rank == 0 && i == 24 && j == 12 && k == 0){
                    Cmpnts ucat_node = corner_arr[k][j][i];
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                }
                */
              //  LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,48,"[Rank %d] Node(k,j,i)=%d,%d,%d finished loops and write.\n", rank, k, j, i);
            }
        }
    }
    PROFILE_FUNCTION_END;
    return 0;
}

InterpolateFieldFromCenterToCorner_Scalar()

PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar	(	PetscReal ***	centfield_arr,
		PetscReal ***	corner_arr,
		UserCtx *	user
	)

Interpolates a scalar field from cell centers to corner nodes.

This version is adapted to write directly into a ghosted local array obtained from DMDAVecGetArray(), which allows using GLOBAL indices for writing to the OWNED portion of the array.

Parameters

[in]	centfield_arr	Input: 3D array (ghosted) of scalar data at cell centers, accessed via GLOBAL indices.
[out]	corner_arr	Output: 3D array (ghosted) where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Parameters

[in]	centfield_arr	Input: 3D array (ghosted) of cell-centered data, accessed via GLOBAL indices.
[out]	corner_arr	Output: 3D array where interpolated node values are stored, also accessed via GLOBAL indices for the owned part.
[in]	user	User context containing DMDA information.

Returns

PetscErrorCode 0 on success.

Definition at line 392 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PetscMPIInt    rank;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xm_node = info.xm, xe_node = xs_node + xm_node;
    PetscInt ys_node = info.ys, ym_node = info.ym, ye_node = ys_node + ym_node;
    PetscInt zs_node = info.zs, zm_node = info.zm, ze_node = zs_node + zm_node;
 
    PetscInt nCellsX = info.mx - 2; // Number of cells in x-direction
    PetscInt nCellsY = info.my - 2; // Number of cells in y-direction
    PetscInt nCellsZ = info.mz - 2; // Number of cells in z-direction
 
 
    // Global grid dimensions (used for valid cell check)
    PetscInt IM = info.mx - 1; // Total nodes in i-direction
    PetscInt JM = info.my - 1; // Total nodes in j-direction
    PetscInt KM = info.mz - 1; // Total nodes in k-direction
 
    // Ghosted array starting indices (from the same info struct)
    PetscInt gxs = info.gxs;
    PetscInt gys = info.gys;
    PetscInt gzs = info.gzs;
 
    PetscInt xs = info.xs;
    PetscInt ys = info.ys;
    PetscInt zs = info.zs;
 
    LOG_ALLOW_SYNC(GLOBAL,LOG_VERBOSE,
        "[Rank %d] Starting -- Node ownership k=%d..%d, j=%d..%d, i=%d..%d\n",
        rank, zs_node, ze_node-1, ys_node, ye_node-1, xs_node, xe_node-1);
 
    // Loop over the GLOBAL indices of the NODES owned by this processor
    for (PetscInt k = zs_node; k < ze_node; k++) {
        for (PetscInt j = ys_node; j < ye_node; j++) {
            for (PetscInt i = xs_node; i < xe_node; i++) {
                PetscReal sum = 0.0;
                PetscInt count = 0;
 
                // DEBUG 1 TEST
                /*
                if(rank == 1 && i == 24 && j == 12 && k == 49){
                    PetscInt i_cell_A = i - 1; 
                    PetscInt i_cell_B = i;
 
                    Cmpnts ucat_A = centfield_arr[k][j][i_cell_A]; // 23
                    Cmpnts ucat_B = centfield_arr[k][j][i_cell_B]; // 24 (out-of-bounds if IM=25)
 
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Read { Valid Cell }Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f) and {Out-of-Bounds Cell}Ucat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i_cell_A, ucat_A.x, ucat_A.y, ucat_A.z,
                             k, j, i_cell_B, ucat_B.x, ucat_B.y, ucat_B.z);
 
                }
                */
                
                // Skip processing the unused last node in each dimension.
                if(i >= IM || j >= JM || k >= KM){
                    continue;
                }
                // Loop over the 8 potential cells surrounding node N(k,j,i) and accumulate values.
                // The index offsets correspond to the relative position of the indices(shifted) that represent cell-centered field values of cells that share the node N(k,j,i) as a corner
                for (PetscInt dk_offset = -1; dk_offset <= 0; dk_offset++) {
                    for (PetscInt dj_offset = -1; dj_offset <= 0; dj_offset++) {
                        for (PetscInt di_offset = -1; di_offset <= 0; di_offset++) {
 
                            // These are still GLOBAL cell indices
                            PetscInt global_cell_k = k + dk_offset;
                            PetscInt global_cell_j = j + dj_offset;
                            PetscInt global_cell_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index
                            if (global_cell_i >= 0 && global_cell_i < nCellsX &&
                                global_cell_j >= 0 && global_cell_j < nCellsY &&
                                global_cell_k >= 0 && global_cell_k < nCellsZ)
                            {
                                    PetscReal cell_val = centfield_arr[global_cell_k + 1][global_cell_j + 1][global_cell_i + 1];
 
                                    LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,49,"[Rank %d] successful read from [%d][%d][%d] -> (%.2f)\n",
                                        rank,global_cell_k,global_cell_j,global_cell_i,cell_val);
                    
                                    sum += cell_val;
                                    count++;                                  
                            }
                        }
                    }
                }
 
                PetscInt i_global_write = i; // Global index in GLOBAL array.
                PetscInt j_global_write = j;
                PetscInt k_global_write = k;
                
                // We write directly into the array using the global loop indices.
                if (count > 0) {
                    corner_arr[k_global_write][j_global_write][i_global_write] = sum / (PetscReal)count;
                } else {
                    // This case should ideally not happen for a valid owned node, but as a failsafe:
                    corner_arr[k_global_write][j_global_write][i_global_write] = 0.0;
                }
               
                // DEBUG 2
                /*
                if(rank == 1){
                    if(i == 11 && j == 11 && k == 49){
                        Cmpnts ucat_node = corner_arr[k][j][i];
                        PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                                 rank, k, j, i,
                                 k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                    }
                }
 
                if(rank == 0 && i == 24 && j == 12 && k == 0){
                    Cmpnts ucat_node = corner_arr[k][j][i];
                    PetscPrintf(PETSC_COMM_WORLD,"[Rank %d] DEBUG TEST at Node(k,j,i)=%d,%d,%d: Wrote CornerUcat[%d][%d][%d]=(%.2f,%.2f,%.2f)\n",
                             rank, k, j, i,
                             k, j, i, ucat_node.x, ucat_node.y, ucat_node.z);
                }
                */
              //  LOG_LOOP_ALLOW_EXACT(LOCAL, LOG_VERBOSE,k,48,"[Rank %d] Node(k,j,i)=%d,%d,%d finished loops and write.\n", rank, k, j, i);
            }
        }
    }
    PROFILE_FUNCTION_END;
    return 0;
}

InterpolateCornerToFaceCenter_Scalar()

PetscErrorCode InterpolateCornerToFaceCenter_Scalar	(	PetscReal ***	corner_arr,
		PetscReal ***	faceX_arr,
		PetscReal ***	faceY_arr,
		PetscReal ***	faceZ_arr,
		UserCtx *	user
	)

Interpolates a scalar field from corner nodes to all face centers.

This routine computes the average of the four corner-node values defining each face of a hexahedral cell:

• X-faces (perpendicular to X): face between (i-1,i) in X-dir
• Y-faces (perpendicular to Y): face between (j-1,j) in Y-dir
• Z-faces (perpendicular to Z): face between (k-1,k) in Z-dir

Parameters

[in]	corner_arr	Ghosted node-centered array (global indexing) from user->fda.
[out]	faceX_arr	Local array for X-faces sized [zm][ym][xm+1].
[out]	faceY_arr	Local array for Y-faces sized [zm][ym+1][xm].
[out]	faceZ_arr	Local array for Z-faces sized [zm+1][ym][xm].
[in]	user	User context containing DMDA 'fda' and GetOwnedCellRange.

Returns

PetscErrorCode 0 on success, non-zero on failure.

Definition at line 2228 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Determine owned-cell ranges based on corner-node ownership
    PetscInt xs, xm, ys, ym, zs, zm;
    ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
 
    // Global exclusive end indices for cells
    PetscInt xe = xs + xm;
    PetscInt ye = ys + ym;
    PetscInt ze = zs + zm;
 
    // --- X‐faces: loops k=zs..ze-1, j=ys..ye-1, i=xs..xe (xm+1 faces per row) ---
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i <= xe; ++i) {
                PetscInt i_loc = i - xs; // 0..xm
                // Average the four corners of the Y-Z face at X = i
                PetscReal sum = corner_arr[k  ][j  ][i]
                              + corner_arr[k+1][j  ][i]
                              + corner_arr[k  ][j+1][i]
                              + corner_arr[k+1][j+1][i];
                faceX_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    // --- Y‐faces: loops k=zs..ze-1, j=ys..ye (ym+1 faces), i=xs..xe-1 ---
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j <= ye; ++j) {
            PetscInt j_loc = j - ys; // 0..ym
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                // Average the four corners of the X-Z face at Y = j
                PetscReal sum = corner_arr[k  ][j][i  ]
                              + corner_arr[k+1][j][i  ]
                              + corner_arr[k  ][j][i+1]
                              + corner_arr[k+1][j][i+1];
                faceY_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    // --- Z‐faces: loops k=zs..ze (zm+1), j=ys..ye-1, i=xs..xe-1 ---
    for (PetscInt k = zs; k <= ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                // Average the four corners of the X-Y face at Z = k
                PetscReal sum = corner_arr[k][j  ][i  ]
                              + corner_arr[k][j  ][i+1]
                              + corner_arr[k][j+1][i  ]
                              + corner_arr[k][j+1][i+1];
                faceZ_arr[k_loc][j_loc][i_loc] = sum * 0.25;
            }
        }
    }
 
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

InterpolateCornerToFaceCenter_Vector()

PetscErrorCode InterpolateCornerToFaceCenter_Vector	(	Cmpnts ***	corner_arr,
		Cmpnts ***	faceX_arr,
		Cmpnts ***	faceY_arr,
		Cmpnts ***	faceZ_arr,
		UserCtx *	user
	)

Interpolates a vector field from corner nodes to all face centers.

Identical to the scalar version, except it averages each component of the Cmpnts struct at the four corner-nodes per face.

Parameters

[in]	corner_arr	Ghosted 3-component array (global node indices).
[out]	faceX_arr	Local array of Cmpnts for X-faces sized [zm][ym][xm+1].
[out]	faceY_arr	Local array of Cmpnts for Y-faces sized [zm][ym+1][xm].
[out]	faceZ_arr	Local array of Cmpnts for Z-faces sized [zm+1][ym][xm].
[in]	user	User context containing DMDA 'fda'.

Returns

PetscErrorCode 0 on success.

Definition at line 2328 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo info;
    PetscMPIInt rank;
 
    PetscFunctionBeginUser;
 
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d starting InterpolateFieldFromCornerToFaceCenter_Vector.\n", rank);
    
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    PetscInt xs, xm, ys, ym, zs, zm;
    ierr = GetOwnedCellRange(&info, 0, &xs, &xm); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 1, &ys, &ym); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info, 2, &zs, &zm); CHKERRQ(ierr);
 
    PetscInt xe = xs + xm;
    PetscInt ye = ys + ym;
    PetscInt ze = zs + zm;
 
    // X-faces
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i <= xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k  ][j  ][i].x + corner_arr[k+1][j  ][i].x
                      + corner_arr[k  ][j+1][i].x + corner_arr[k+1][j+1][i].x;
                sum.y = corner_arr[k  ][j  ][i].y + corner_arr[k+1][j  ][i].y
                      + corner_arr[k  ][j+1][i].y + corner_arr[k+1][j+1][i].y;
                sum.z = corner_arr[k  ][j  ][i].z + corner_arr[k+1][j  ][i].z
                      + corner_arr[k  ][j+1][i].z + corner_arr[k+1][j+1][i].z;
                faceX_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceX_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceX_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d x-face Interpolation complete.\n", rank);
    
    // Y-faces
    for (PetscInt k = zs; k < ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j <= ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k  ][j][i  ].x + corner_arr[k+1][j][i  ].x
                      + corner_arr[k  ][j][i+1].x + corner_arr[k+1][j][i+1].x;
                sum.y = corner_arr[k  ][j][i  ].y + corner_arr[k+1][j][i  ].y
                      + corner_arr[k  ][j][i+1].y + corner_arr[k+1][j][i+1].y;
                sum.z = corner_arr[k  ][j][i  ].z + corner_arr[k+1][j][i  ].z
                      + corner_arr[k  ][j][i+1].z + corner_arr[k+1][j][i+1].z;
                faceY_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceY_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceY_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d y-face Interpolation complete.\n", rank);
    
    // Z-faces
    for (PetscInt k = zs; k <= ze; ++k) {
        PetscInt k_loc = k - zs;
        for (PetscInt j = ys; j < ye; ++j) {
            PetscInt j_loc = j - ys;
            for (PetscInt i = xs; i < xe; ++i) {
                PetscInt i_loc = i - xs;
                Cmpnts sum = {0,0,0};
                sum.x = corner_arr[k][j  ][i  ].x + corner_arr[k][j  ][i+1].x
                      + corner_arr[k][j+1][i  ].x + corner_arr[k][j+1][i+1].x;
                sum.y = corner_arr[k][j  ][i  ].y + corner_arr[k][j  ][i+1].y
                      + corner_arr[k][j+1][i  ].y + corner_arr[k][j+1][i+1].y;
                sum.z = corner_arr[k][j  ][i  ].z + corner_arr[k][j  ][i+1].z
                      + corner_arr[k][j+1][i  ].z + corner_arr[k][j+1][i+1].z;
                faceZ_arr[k_loc][j_loc][i_loc].x = sum.x * 0.25;
                faceZ_arr[k_loc][j_loc][i_loc].y = sum.y * 0.25;
                faceZ_arr[k_loc][j_loc][i_loc].z = sum.z * 0.25;
            }
        }
    }
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_TRACE,
                   "Rank %d z-face Interpolation complete.\n", rank);
    
    PROFILE_FUNCTION_END;
    PetscFunctionReturn(0);
}

Here is the call graph for this function: