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	InterpolateFieldFromCenterToCorner_Scalar (PetscReal ***field_arr, PetscReal ***centfield_arr, UserCtx *user)
	Interpolates a scalar field from cell centers to corner nodes.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Vector (Cmpnts ***field_arr, Cmpnts ***centfield_arr, UserCtx *user)
	Interpolates a vector field from cell centers to corner nodes.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Vector_Petsc (Cmpnts ***centfield_arr, Cmpnts ***corner_arr, UserCtx *user)
	Interpolates a vector field from cell centers to corner nodes.
PetscErrorCode	InterpolateFieldFromCenterToCorner_Scalar_Petsc (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	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_Petsc((PetscReal***)(centfield_ptr), (PetscReal***)(corner_ptr), (user_ctx)) : \
        InterpolateFieldFromCenterToCorner_Vector_Petsc((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 64 of file interpolation.h.

                                                                                                                               : \
        InterpolateFieldFromCenterToCorner_Vector_Petsc((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 87 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 109 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 130 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 622 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 702 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.

Definition at line 958 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;
    PetscInt  *pids = NULL; // For logging particle IDs in warnings
    void      *swarmOut = 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);
 
    /* STAGE 1: Center-to-Corner Calculation with Communication */
 
    // (A) Create the corner vectors.
    if(user->CellFieldAtCorner == NULL){
        LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Creating new global corner vector for field '%s'.\n", rank, fieldName);
        ierr = DMCreateGlobalVector(fda, &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(fda, &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(fda, 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(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);
    
    size = 0;
 
    ierr = VecGetLocalSize(cornerGlobal, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Corner global vector LOCAL size for field '%s': %d.\n", rank, fieldName, (PetscInt)size);
    
    size = 0;
 
    ierr = VecGetLocalSize(fieldLocal_cellCentered, &size); CHKERRQ(ierr);
    LOG_ALLOW(LOCAL, LOG_TRACE, "[Rank %d] Cell-centered local vector LOCAL 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);
 
    PetscInt xs,ys,zs,gxs,gys,gzs;
    ierr  = DMDAGetCorners(fda,&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(fda,&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(fda, 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(fda, cornerGlobal, &cornerPtr_write); CHKERRQ(ierr);
    
    ierr = DMDAVecRestoreArrayRead(fda, 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(fda, cornerGlobal, INSERT_VALUES, cornerLocal); CHKERRQ(ierr);
    ierr = DMGlobalToLocalEnd(fda, 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 FDA on all ranks.\n");
       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(fda, 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, 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++) {
        PetscInt iCell_global = cellIDs[3 * p + 0];
        PetscInt jCell_global = cellIDs[3 * p + 1];
        PetscInt kCell_global = cellIDs[3 * p + 2];
 
 
        LOG_ALLOW(LOCAL, LOG_VERBOSE,
          "[Rank %d] Particle PID %lld: global cell=(%d,%d,%d), weights=(%.4f,%.4f,%.4f)\n",
          rank, (long long)pids[p], iCell_global, jCell_global, kCell_global,
          weights[3 * p + 0], weights[3 * p + 1], weights[3 * p + 2]);
        // Safety check: Ensure the entire 8-node stencil is within the valid memory
        // region (owned + ghosts) of the local array.
 
        if (iCell_global < info.gxs || iCell_global >= info.gxs + info.gxm - 1 ||
            jCell_global < info.gxs || jCell_global >= info.gys + info.gym - 1 ||
            kCell_global < info.gxs || kCell_global >= 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)pids[p], iCell_global, jCell_global, kCell_global, 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;
        }
 
        PetscReal a1 = weights[3 * p + 0];
        PetscReal a2 = weights[3 * p + 1];
        PetscReal a3 = weights[3 * p + 2];
 
        ierr = InterpolateEulerFieldToSwarmForParticle(
                  fieldName,
                  cornerPtr_read_with_ghosts,
                  iCell_global, jCell_global, kCell_global,
                  a1, a2, a3,
                  swarmOut, p, bs);
        CHKERRQ(ierr);
    }
 
    /* (H) Restore all retrieved arrays and vectors */
    ierr = DMDAVecRestoreArrayRead(fda, cornerLocal, &cornerPtr_read_with_ghosts); CHKERRQ(ierr);
    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);
 
    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 1276 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 *GLOBAL* Vector to the function below! 
  ierr = InterpolateEulerFieldToSwarm(user, user->lUcat, 
                                      "Ucat", 
                                      "velocity"); CHKERRQ(ierr);
  /* 
     2) (OPTIONAL) If you have more fields, you can Interpolatete them similarly.
 
     For example, if user->Tcat is a scalar Vec for "temperature" and the DMSwarm
     has a field "Temperature" (with dof=1), you could do:
 
         ierr = InterpolateOneFieldOverSwarm(user, user->Tcat,
                                             "temperature",
                                             "swarmTemperature"); CHKERRQ(ierr);
 
     For pressure:
         ierr = InterpolateOneFieldOverSwarm(user, user->Pcat,
                                             "pressure",
                                             "swarmPressure"); CHKERRQ(ierr);
     
     ...and so forth.
   */
 
  /* 
     3) Optionally, synchronize or log that all fields are done 
  */
  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.

This function estimates the value of a scalar field at the center of each computational cell (whose origin node is owned by this rank) by averaging the values from the 8 corner nodes defining that cell.

Parameters

[in]	field_arr	Input: 3D array view (ghosted) of scalar data at nodes, obtained from a Vec associated with user->da. Accessed via GLOBAL node indices. Ghost values must be up-to-date.
[out]	centfield_arr	Output: 3D local array (0-indexed) where interpolated cell center values are stored. Sized [zm_cell_local][ym_cell_local][xm_cell_local] by the caller. centfield_arr[k_loc][j_loc][i_loc] stores the value for the cell corresponding to the (i_loc, j_loc, k_loc)-th owned cell origin.
[in]	user	User context containing DMDA information (da, IM, JM, KM).

Returns

PetscErrorCode 0 on success.

Definition at line 166 of file interpolation.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
                   "Rank %d starting InterpolateFieldFromCornerToCenter_Scalar.\n", rank);
 
    // Get local info based on user->da (which is node-based but defines cell origins)
    DMDALocalInfo info_da_nodes;
    ierr = DMDAGetLocalInfo(user->da, &info_da_nodes); CHKERRQ(ierr);
 
    // Determine owned CELL ranges based on owning their origin NODES (from user->da)
    PetscInt xs_cell_global_i, xm_cell_local_i; // Global start cell index, local count of owned cells
    PetscInt ys_cell_global_j, ym_cell_local_j;
    PetscInt zs_cell_global_k, zm_cell_local_k;
 
    ierr = GetOwnedCellRange(&info_da_nodes, 0, &xs_cell_global_i, &xm_cell_local_i); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info_da_nodes, 1, &ys_cell_global_j, &ym_cell_local_j); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info_da_nodes, 2, &zs_cell_global_k, &zm_cell_local_k); CHKERRQ(ierr);
 
    // Exclusive end global cell indices
    PetscInt xe_cell_global_i_excl = xs_cell_global_i + xm_cell_local_i;
    PetscInt ye_cell_global_j_excl = ys_cell_global_j + ym_cell_local_j;
    PetscInt ze_cell_global_k_excl = zs_cell_global_k + zm_cell_local_k;
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d: Processing Owned Cells (global indices) k=%d..%d (count %d), j=%d..%d (count %d), i=%d..%d (count %d)\n",
              rank,
              zs_cell_global_k, ze_cell_global_k_excl -1, zm_cell_local_k,
              ys_cell_global_j, ye_cell_global_j_excl -1, ym_cell_local_j,
              xs_cell_global_i, xe_cell_global_i_excl -1, xm_cell_local_i);
 
    // Loop over the GLOBAL indices of the CELLS whose origin nodes are owned by this processor
    // (according to user->da partitioning)
    for (PetscInt k_glob_cell = zs_cell_global_k; k_glob_cell < ze_cell_global_k_excl; k_glob_cell++) {
        PetscInt k_local = k_glob_cell - zs_cell_global_k; // 0-based local index for centfield_arr
 
        for (PetscInt j_glob_cell = ys_cell_global_j; j_glob_cell < ye_cell_global_j_excl; j_glob_cell++) {
            PetscInt j_local = j_glob_cell - ys_cell_global_j; // 0-based local index
 
            for (PetscInt i_glob_cell = xs_cell_global_i; i_glob_cell < xe_cell_global_i_excl; i_glob_cell++) {
                PetscInt i_local = i_glob_cell - xs_cell_global_i; // 0-based local index
 
                PetscReal sum = 0.0;
                // Count is always 8 for a hexahedral cell if all nodes are accessible
                // PetscInt count = 0;
 
                // Loop over the 8 NODE indices that define cell C(i_glob_cell, j_glob_cell, k_glob_cell)
                for (PetscInt dk = 0; dk < 2; dk++) {
                    for (PetscInt dj = 0; dj < 2; dj++) {
                        for (PetscInt di = 0; di < 2; di++) {
                            PetscInt ni_glob = i_glob_cell + di; // Global NODE index for corner
                            PetscInt nj_glob = j_glob_cell + dj; // Global NODE index for corner
                            PetscInt nk_glob = k_glob_cell + dk; // Global NODE index for corner
 
                            // Access the input 'field_arr' (node values) using GLOBAL node indices.
                            sum += field_arr[nk_glob][nj_glob][ni_glob];
                            // count++;
                        }
                    }
                }
 
                PetscReal center_value = sum / 8.0;
 
                // Store the calculated cell center value into the output array 'centfield_arr'
                // using the LOCAL 0-based cell index.
                // Defensive check (should not be strictly necessary if caller allocated centfield_arr correctly)
                if (i_local < 0 || i_local >= xm_cell_local_i ||
                    j_local < 0 || j_local >= ym_cell_local_j ||
                    k_local < 0 || k_local >= zm_cell_local_k) {
                    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local index out of bounds for centfield_arr write in Scalar version!");
                }
                centfield_arr[k_local][j_local][i_local] = center_value;
 
            } // End loop i_glob_cell
        } // End loop j_glob_cell
    } // End loop k_glob_cell
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
                   "Rank %d completed InterpolateFieldFromCornerToCenter_Scalar.\n", rank);
    
    PROFILE_FUNCTION_END;
              
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

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.

This function estimates the value of a vector field at the center of each computational cell (whose origin node is owned by this rank) by averaging the values from the 8 corner nodes defining that cell.

Parameters

[in]	field_arr	Input: 3D array view (ghosted) of vector data at nodes, obtained from a Vec associated with user->fda. Accessed via GLOBAL node indices. Ghost values must be up-to-date.
[out]	centfield_arr	Output: 3D local array (0-indexed) where interpolated cell center values are stored. Sized [zm_cell_local][ym_cell_local][xm_cell_local] by the caller. centfield_arr[k_loc][j_loc][i_loc] stores the value for the cell corresponding to the (i_loc, j_loc, k_loc)-th owned cell origin.
[in]	user	User context containing DMDA information (fda, IM, JM, KM).

Returns

PetscErrorCode 0 on success.

Definition at line 39 of file interpolation.c.

{
    PetscErrorCode ierr;
    PetscMPIInt    rank;
    PetscFunctionBeginUser;
    PROFILE_FUNCTION_BEGIN;
 
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
                   "Rank %d starting InterpolateFieldFromCornerToCenter_Vector.\n", rank);
 
    // Get local info based on the NODE-based DMDA (user->fda)
    // This DM defines the ownership of the nodes from which we interpolate,
    // and thus defines the cells whose centers we will compute.
    DMDALocalInfo info_fda;
    ierr = DMDAGetLocalInfo(user->fda, &info_fda); CHKERRQ(ierr);
 
    // Determine owned CELL ranges based on owning their origin NODES (from user->fda)
    PetscInt xs_cell_global_i, xm_cell_local_i; // Global start cell index, local count of owned cells
    PetscInt ys_cell_global_j, ym_cell_local_j;
    PetscInt zs_cell_global_k, zm_cell_local_k;
 
    
    ierr = GetOwnedCellRange(&info_fda, 0, &xs_cell_global_i, &xm_cell_local_i); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info_fda, 1, &ys_cell_global_j, &ym_cell_local_j); CHKERRQ(ierr);
    ierr = GetOwnedCellRange(&info_fda, 2, &zs_cell_global_k, &zm_cell_local_k); CHKERRQ(ierr);
 
    // Exclusive end global cell indices
    PetscInt xe_cell_global_i_excl = xs_cell_global_i + xm_cell_local_i;
    PetscInt ye_cell_global_j_excl = ys_cell_global_j + ym_cell_local_j;
    PetscInt ze_cell_global_k_excl = zs_cell_global_k + zm_cell_local_k;
 
    LOG_ALLOW(LOCAL, LOG_TRACE, "Rank %d : Processing Owned Cells (global indices) k=%d..%d (count %d), j=%d..%d (count %d), i=%d..%d (count %d)\n",
              rank,
              zs_cell_global_k, ze_cell_global_k_excl -1, zm_cell_local_k,
              ys_cell_global_j, ye_cell_global_j_excl -1, ym_cell_local_j,
              xs_cell_global_i, xe_cell_global_i_excl -1, xm_cell_local_i);
 
    // Loop over the GLOBAL indices of the CELLS whose origin nodes are owned by this processor
    // (according to user->fda partitioning)
    for (PetscInt k_glob_cell = zs_cell_global_k; k_glob_cell < ze_cell_global_k_excl; k_glob_cell++) {
        PetscInt k_local = k_glob_cell - zs_cell_global_k; // 0-based local index for centfield_arr
 
        for (PetscInt j_glob_cell = ys_cell_global_j; j_glob_cell < ye_cell_global_j_excl; j_glob_cell++) {
            PetscInt j_local = j_glob_cell - ys_cell_global_j; // 0-based local index
 
            for (PetscInt i_glob_cell = xs_cell_global_i; i_glob_cell < xe_cell_global_i_excl; i_glob_cell++) {
                PetscInt i_local = i_glob_cell - xs_cell_global_i; // 0-based local index
 
                Cmpnts sum = {0.0, 0.0, 0.0};
                // Count is always 8 for a hexahedral cell
                // PetscInt count = 0; // Not strictly needed if we assume 8 corners
 
                // Loop over the 8 corner NODES that define cell C(i_glob_cell, j_glob_cell, k_glob_cell)
                for (PetscInt dk = 0; dk < 2; dk++) {
                    for (PetscInt dj = 0; dj < 2; dj++) {
                        for (PetscInt di = 0; di < 2; di++) {
                            PetscInt ni_glob = i_glob_cell + di; // Global NODE index i or i+1
                            PetscInt nj_glob = j_glob_cell + dj; // Global NODE index j or j+1
                            PetscInt nk_glob = k_glob_cell + dk; // Global NODE index k or k+1
 
                            // Access the input 'field_arr' (node values) using GLOBAL node indices.
                            // Ghosts in field_arr must be valid.
                            // DMDAVecGetArrayRead handles mapping global indices to local memory.
                            sum.x += field_arr[nk_glob][nj_glob][ni_glob].x;
                            sum.y += field_arr[nk_glob][nj_glob][ni_glob].y;
                            sum.z += field_arr[nk_glob][nj_glob][ni_glob].z;
                            // count++;
                        }
                    }
                }
 
                // Calculate average value representing the cell center
                Cmpnts center_value;
                center_value.x = sum.x / 8.0;
                center_value.y = sum.y / 8.0;
                center_value.z = sum.z / 8.0;
 
                // Store result into the local, 0-indexed output array 'centfield_arr'
                // Defensive check (should not be strictly necessary if caller allocated centfield_arr
                // with dimensions xm_cell_local_i, ym_cell_local_j, zm_cell_local_k)
                if (i_local < 0 || i_local >= xm_cell_local_i ||
                    j_local < 0 || j_local >= ym_cell_local_j ||
                    k_local < 0 || k_local >= zm_cell_local_k) {
                    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Local index out of bounds for centfield_arr write!");
                }
                centfield_arr[k_local][j_local][i_local] = center_value;
 
            } // End loop i_glob_cell
        } // End loop j_glob_cell
    } // End loop k_glob_cell
 
    LOG_ALLOW_SYNC(GLOBAL, LOG_DEBUG,
                   "Rank %d completed InterpolateFieldFromCornerToCenter_Vector.\n", rank);
    
    PROFILE_FUNCTION_END;
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

InterpolateFieldFromCenterToCorner_Scalar()

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

Interpolates a scalar field from cell centers to corner nodes.

This function estimates the value of a scalar field at each grid node by averaging the values from the cell centers of the cells surrounding that node (up to 8). It handles physical boundaries by averaging only the available adjacent cells.

Assumes input centfield_arr is from a ghosted local vector associated with user->da (DOF=1, s=2) and output field_arr is from a ghosted local vector also associated with user->da (DOF=1, s=2). Input array uses GLOBAL cell indices, output array uses GLOBAL node indices.

Parameters

[in]	centfield_arr	Input: 3D array (ghosted) of scalar data at cell centers, accessed via GLOBAL cell indices (k=0..KM-1, j=0..JM-1, i=0..IM-1).
[out]	field_arr	Output: 3D array (ghosted) where interpolated node values are stored, accessed via GLOBAL node indices (k=0..KM, j=0..JM, i=0..IM).
[in]	user	User context containing DMDA information (da).

Returns

PetscErrorCode 0 on success.

InterpolateFieldFromCenterToCorner_Vector()

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

Interpolates a vector field from cell centers to corner nodes.

This function estimates the value of a vector field at each grid node by averaging the vector values from the cell centers of the cells surrounding that node (up to 8). It handles physical boundaries by averaging only the available adjacent cells.

Assumes input centfield_arr is from a ghosted local vector (e.g., representing ucat, stored using node-indexing convention) and output field_arr is a ghosted local vector associated with user->fda (DOF=3, s=2), accessed using global node indices.

Parameters

[in]	centfield_arr	Input: 3D array (ghosted) of vector data conceptually at cell centers, accessed via GLOBAL indices respecting the storage convention (e.g., ucat[k][j][i] uses node index i but represents cell C(i,j,k) for interior).
[out]	field_arr	Output: 3D array (ghosted) where interpolated node values are stored, accessed via GLOBAL node indices (k=0..KM, j=0..JM, i=0..IM).
[in]	user	User context containing DMDA information (da and fda).

Returns

PetscErrorCode 0 on success.

InterpolateFieldFromCenterToCorner_Vector_Petsc()

PetscErrorCode InterpolateFieldFromCenterToCorner_Vector_Petsc	(	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 270 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 InterpolateFieldFromCenterToCorner_Vector_Petsc: 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_Petsc()

PetscErrorCode InterpolateFieldFromCenterToCorner_Scalar_Petsc	(	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.

Definition at line 424 of file interpolation.c.

{
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    PROFILE_FUNCTION_BEGIN;
    ierr = DMDAGetLocalInfo(user->fda, &info); CHKERRQ(ierr);
 
    // Node ownership range (GLOBAL indices)
    PetscInt xs_node = info.xs, xe_node = info.xs + info.xm;
    PetscInt ys_node = info.ys, ye_node = info.ys + info.ym;
    PetscInt zs_node = info.zs, ze_node = info.zs + info.zm;
 
    // Global grid dimensions (number of nodes)
    PetscInt IM = info.mx-1;
    PetscInt JM = info.my-1;
    PetscInt KM = info.mz-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;
 
                // Loop over the 8 potential cells surrounding node N(k,j,i)
                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++) {
                            // This is the global index of the cell-centered data point
                            PetscInt cell_idx_k = k + dk_offset;
                            PetscInt cell_idx_j = j + dj_offset;
                            PetscInt cell_idx_i = i + di_offset;
 
                            // Check if this corresponds to a valid GLOBAL cell index.
                            // Cell indices run from 0 to NumberOfNodes-2.
                            if (cell_idx_i >= 0 && cell_idx_i < IM - 1 &&
                                cell_idx_j >= 0 && cell_idx_j < JM - 1 &&
                                cell_idx_k >= 0 && cell_idx_k < KM - 1)
                            {
                                sum += centfield_arr[cell_idx_k][cell_idx_j][cell_idx_i];
                                count++;
                            }
                        }
                    }
                }
 
                // Calculate average and write to the output array using GLOBAL indices
                if (count > 0) {
                    corner_arr[k][j][i] = sum / (PetscReal)count;
                } else {
                    corner_arr[k][j][i] = 0.0; // Failsafe
                }
            }
        }
    }
    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 2004 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 2104 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: