Field3D
MACFieldUtil.h File Reference

Contains utility functions for MAC fields, such as conversion to cell-centered fields. More...

#include "MACField.h"
#include "ns.h"

Go to the source code of this file.

Functions

template<class Field_T , class Data_T >
void convertCellCenteredToMAC (typename Field_T::Ptr cc, typename MACField< Data_T >::Ptr mac)
 Converts the cell-centered field to a MAC field.
template<class Data_T , class Field_T >
FIELD3D_NAMESPACE_OPEN void convertMACToCellCentered (typename MACField< Data_T >::Ptr mac, typename Field_T::Ptr cc)
 Converts the MAC field to a cell-centered field.

Detailed Description

Contains utility functions for MAC fields, such as conversion to cell-centered fields.

Definition in file MACFieldUtil.h.


Function Documentation

template<class Data_T , class Field_T >
void convertMACToCellCentered ( typename MACField< Data_T >::Ptr  mac,
typename Field_T::Ptr  cc 
)

Converts the MAC field to a cell-centered field.

Sets up the cell-centered target field given a MAC field.

Definition at line 78 of file MACFieldUtil.h.

References FieldRes::dataWindow(), FieldRes::extents(), field_dynamic_cast(), MatrixFieldMapping::localToWorld(), FieldRes::mapping(), and MACField< Data_T >::value().

{
  // Make sure the extents and data window match
  if (cc->extents().min != mac->extents().min ||
      cc->extents().max != mac->extents().max ||
      cc->dataWindow().min != mac->dataWindow().min ||
      cc->dataWindow().max != mac->dataWindow().max ) {
    cc->setSize(mac->extents(), mac->dataWindow());
  }

  // Make sure mapping matches
  if (!cc->mapping()->isIdentical(mac->mapping())) {
    cc->setMapping(mac->mapping());
  }

  // MAC velocities are in simulation space (axis-aligned to the
  // mapping) because the values are stored on the faces, so rotate
  // vectors from simulation-space to world-space when transferring
  // from MAC to cell-centered

  bool rotateVector = false;
  M44d ssToWsMtx;
  MatrixFieldMapping::Ptr mapping =
    field_dynamic_cast<MatrixFieldMapping>(mac->mapping());
  if (mapping) {
    M44d localToWorldMtx = mapping->localToWorld();
    V3d scale, rot, trans, shear;
    if (extractSHRT(localToWorldMtx, scale, shear, rot, trans, false)) {
      ssToWsMtx.rotate(rot);
      if (rot.length2() > FLT_EPSILON)
        rotateVector = true;
    }
  }

  // Loop over all the voxels in the output field ---

  typename Field_T::iterator i = cc->begin();
  typename Field_T::iterator end = cc->end();

  if (rotateVector) {
    for (; i != end; ++i) {
      *i = mac->value(i.x, i.y, i.z) * ssToWsMtx;
    }
  } else {
    for (; i != end; ++i) {
      *i = mac->value(i.x, i.y, i.z);
    }
  }
}
template<class Field_T , class Data_T >
void convertCellCenteredToMAC ( typename Field_T::Ptr  cc,
typename MACField< Data_T >::Ptr  mac 
)

Converts the cell-centered field to a MAC field.

Definition at line 132 of file MACFieldUtil.h.

References FieldRes::dataWindow(), FieldRes::extents(), FIELD3D_EXTRACT_SHRT, field_dynamic_cast(), MatrixFieldMapping::localToWorld(), FieldRes::mapping(), FieldRes::setMapping(), ResizableField< Data_T >::setSize(), MACField< Data_T >::u(), MACField< Data_T >::v(), and MACField< Data_T >::w().

{
  // Make sure the extents and data window match
  if (mac->extents().min != cc->extents().min ||
      mac->extents().max != cc->extents().max ||
      mac->dataWindow().min != cc->dataWindow().min ||
      mac->dataWindow().max != cc->dataWindow().max ) {
    mac->setSize(cc->extents(), cc->dataWindow());
  }

  // Make sure mapping matches
  if (!mac->mapping()->isIdentical(cc->mapping())) {
    mac->setMapping(cc->mapping());
  }
  
  Box3i data = mac->dataWindow();

  // MAC velocities are in simulation space (axis-aligned to the
  // mapping) because the values are stored on the faces, so rotate
  // vectors from world-space to simulation-space when transferring
  // from cell-centered to MAC

  bool rotateVector = false;
  M44d wsToSsMtx;
  MatrixFieldMapping::Ptr mapping =
    field_dynamic_cast<MatrixFieldMapping>(mac->mapping());
  if (mapping) {
    M44d localToWorld = mapping->localToWorld();
    V3d scale, rot, trans, shear;
    if (FIELD3D_EXTRACT_SHRT(localToWorld, scale, shear, rot, trans, false)) {
      wsToSsMtx.rotate(-rot);
      rotateVector = true;
    }
  }

  // Use a pointer to a field below so we can substitute it out for
  // our intermediate, rotated field if necessary, without needing to
  // duplicate the loops.  This should be more efficient CPU-wise so
  // we don't need to do 3 matrix multiplies per cell-centered voxel
  // because it's used in 3 separate loops (1 per MAC face).
  typename Field_T::Ptr src = cc;

  typename Field_T::Ptr ccRotated;
  if (rotateVector) {
    ccRotated =
      typename Field_T::Ptr(new Field_T);
    ccRotated->matchDefinition(cc);

    typename Field_T::const_iterator iIn = cc->cbegin();
    typename Field_T::const_iterator endIn = cc->cend();
    typename Field_T::iterator iOut = ccRotated->begin();

    for (; iIn != endIn; ++iIn, ++iOut) {
      *iOut = *iIn * wsToSsMtx;
    }
    src = ccRotated;
  }

  // Set the u edge value to their closest voxel
  for (int k = data.min.z; k <= data.max.z; k++) {
    for (int j = data.min.y; j <= data.max.y; j++) {
      mac->u(data.min.x, j, k) = src->value(data.min.x, j, k).x;
      mac->u(data.max.x + 1, j, k) = src->value(data.max.x, j, k).x;
    }
  }

  // Set the v edge value to their closest voxel
  for (int k = data.min.z; k <= data.max.z; k++) {
    for (int i = data.min.x; i <= data.max.x; i++) {
      mac->v(i, data.min.y, k) = src->value(i, data.min.y, k).y;
      mac->v(i, data.max.y + 1, k) = src->value(i, data.max.y, k).y;
    }
  }

  // Set the w edge value to their closest voxel
  for (int j = data.min.y; j <= data.max.y; j++) {
    for (int i = data.min.x; i <= data.max.x; i++) {
      mac->w(i, j, data.min.z) = src->value(i, j, data.min.z).z;
      mac->w(i, j, data.max.z + 1) = src->value(i, j, data.max.z).z;
    }
  }

  // Loop over internal u values
  for (int k = data.min.z; k <= data.max.z; ++k) {
    for (int j = data.min.y; j <= data.max.y; ++j) {
      for (int i = data.min.x + 1; i <= data.max.x; ++i) {
        mac->u(i, j, k) = 
          (src->value(i, j, k).x + src->value(i - 1, j, k).x) * 0.5;
      }
    }
  }

  // Loop over internal v values
  for (int k = data.min.z; k <= data.max.z; ++k) {
    for (int j = data.min.y + 1; j <= data.max.y; ++j) {
      for (int i = data.min.x; i <= data.max.x; ++i) {
        mac->v(i, j, k) = 
          (src->value(i, j, k).y + src->value(i, j - 1, k).y) * 0.5;
      }
    }
  }

  // Loop over internal w values
  for (int k = data.min.z + 1; k <= data.max.z; ++k) {
    for (int j = data.min.y; j <= data.max.y; ++j) {
      for (int i = data.min.x; i <= data.max.x; ++i) {
        mac->w(i, j, k) = 
          (src->value(i, j, k).z + src->value(i, j, k - 1).z) * 0.5;
      }
    }
  }
}