Field3D
MACFieldUtil.h
Go to the documentation of this file.
00001 //----------------------------------------------------------------------------//
00002 
00003 /*
00004  * Copyright (c) 2009 Sony Pictures Imageworks Inc
00005  *
00006  * All rights reserved.
00007  *
00008  * Redistribution and use in source and binary forms, with or without
00009  * modification, are permitted provided that the following conditions
00010  * are met:
00011  *
00012  * Redistributions of source code must retain the above copyright
00013  * notice, this list of conditions and the following disclaimer.
00014  * Redistributions in binary form must reproduce the above copyright
00015  * notice, this list of conditions and the following disclaimer in the
00016  * documentation and/or other materials provided with the
00017  * distribution.  Neither the name of Sony Pictures Imageworks nor the
00018  * names of its contributors may be used to endorse or promote
00019  * products derived from this software without specific prior written
00020  * permission.
00021  *
00022  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00023  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00024  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
00025  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE
00026  * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
00027  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00028  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
00029  * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00030  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
00031  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
00032  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
00033  * OF THE POSSIBILITY OF SUCH DAMAGE.
00034  */
00035 
00036 //----------------------------------------------------------------------------//
00037 
00043 //----------------------------------------------------------------------------//
00044 
00045 #ifndef _INCLUDED_Field3D_MACFieldUtil_H_
00046 #define _INCLUDED_Field3D_MACFieldUtil_H_
00047 
00048 #include "MACField.h"
00049 
00050 //----------------------------------------------------------------------------//
00051 
00052 #include "ns.h"
00053 
00054 FIELD3D_NAMESPACE_OPEN
00055 
00056 //----------------------------------------------------------------------------//
00057 // Utility functions
00058 //----------------------------------------------------------------------------//
00059 
00061 template <class Data_T, class Field_T>
00062 void convertMACToCellCentered(typename MACField<Data_T>::Ptr mac,
00063                               typename Field_T::Ptr cc);
00064 
00065 //----------------------------------------------------------------------------//
00066 
00068 template <class Field_T, class Data_T>
00069 void convertCellCenteredToMAC(typename Field_T::Ptr cc,
00070                               typename MACField<Data_T>::Ptr mac);
00071 
00072 //----------------------------------------------------------------------------//
00073 // Implementations
00074 //----------------------------------------------------------------------------//
00075 
00077 template <class Data_T, class Field_T>
00078 void convertMACToCellCentered(typename MACField<Data_T>::Ptr mac,
00079                               typename Field_T::Ptr cc)
00080 {
00081   // Make sure the extents and data window match
00082   if (cc->extents().min != mac->extents().min ||
00083       cc->extents().max != mac->extents().max ||
00084       cc->dataWindow().min != mac->dataWindow().min ||
00085       cc->dataWindow().max != mac->dataWindow().max ) {
00086     cc->setSize(mac->extents(), mac->dataWindow());
00087   }
00088 
00089   // Make sure mapping matches
00090   if (!cc->mapping()->isIdentical(mac->mapping())) {
00091     cc->setMapping(mac->mapping());
00092   }
00093 
00094   // MAC velocities are in simulation space (axis-aligned to the
00095   // mapping) because the values are stored on the faces, so rotate
00096   // vectors from simulation-space to world-space when transferring
00097   // from MAC to cell-centered
00098 
00099   bool rotateVector = false;
00100   M44d ssToWsMtx;
00101   MatrixFieldMapping::Ptr mapping =
00102     field_dynamic_cast<MatrixFieldMapping>(mac->mapping());
00103   if (mapping) {
00104     M44d localToWorldMtx = mapping->localToWorld();
00105     V3d scale, rot, trans, shear;
00106     if (extractSHRT(localToWorldMtx, scale, shear, rot, trans, false)) {
00107       ssToWsMtx.rotate(rot);
00108       if (rot.length2() > FLT_EPSILON)
00109         rotateVector = true;
00110     }
00111   }
00112 
00113   // Loop over all the voxels in the output field ---
00114 
00115   typename Field_T::iterator i = cc->begin();
00116   typename Field_T::iterator end = cc->end();
00117 
00118   if (rotateVector) {
00119     for (; i != end; ++i) {
00120       *i = mac->value(i.x, i.y, i.z) * ssToWsMtx;
00121     }
00122   } else {
00123     for (; i != end; ++i) {
00124       *i = mac->value(i.x, i.y, i.z);
00125     }
00126   }
00127 }
00128 
00129 //----------------------------------------------------------------------------//
00130 
00131 template <class Field_T, class Data_T>
00132 void convertCellCenteredToMAC(typename Field_T::Ptr cc,
00133                               typename MACField<Data_T>::Ptr mac)
00134 {
00135   // Make sure the extents and data window match
00136   if (mac->extents().min != cc->extents().min ||
00137       mac->extents().max != cc->extents().max ||
00138       mac->dataWindow().min != cc->dataWindow().min ||
00139       mac->dataWindow().max != cc->dataWindow().max ) {
00140     mac->setSize(cc->extents(), cc->dataWindow());
00141   }
00142 
00143   // Make sure mapping matches
00144   if (!mac->mapping()->isIdentical(cc->mapping())) {
00145     mac->setMapping(cc->mapping());
00146   }
00147   
00148   Box3i data = mac->dataWindow();
00149 
00150   // MAC velocities are in simulation space (axis-aligned to the
00151   // mapping) because the values are stored on the faces, so rotate
00152   // vectors from world-space to simulation-space when transferring
00153   // from cell-centered to MAC
00154 
00155   bool rotateVector = false;
00156   M44d wsToSsMtx;
00157   MatrixFieldMapping::Ptr mapping =
00158     field_dynamic_cast<MatrixFieldMapping>(mac->mapping());
00159   if (mapping) {
00160     M44d localToWorld = mapping->localToWorld();
00161     V3d scale, rot, trans, shear;
00162     if (FIELD3D_EXTRACT_SHRT(localToWorld, scale, shear, rot, trans, false)) {
00163       wsToSsMtx.rotate(-rot);
00164       rotateVector = true;
00165     }
00166   }
00167 
00168   // Use a pointer to a field below so we can substitute it out for
00169   // our intermediate, rotated field if necessary, without needing to
00170   // duplicate the loops.  This should be more efficient CPU-wise so
00171   // we don't need to do 3 matrix multiplies per cell-centered voxel
00172   // because it's used in 3 separate loops (1 per MAC face).
00173   typename Field_T::Ptr src = cc;
00174 
00175   typename Field_T::Ptr ccRotated;
00176   if (rotateVector) {
00177     ccRotated =
00178       typename Field_T::Ptr(new Field_T);
00179     ccRotated->matchDefinition(cc);
00180 
00181     typename Field_T::const_iterator iIn = cc->cbegin();
00182     typename Field_T::const_iterator endIn = cc->cend();
00183     typename Field_T::iterator iOut = ccRotated->begin();
00184 
00185     for (; iIn != endIn; ++iIn, ++iOut) {
00186       *iOut = *iIn * wsToSsMtx;
00187     }
00188     src = ccRotated;
00189   }
00190 
00191   // Set the u edge value to their closest voxel
00192   for (int k = data.min.z; k <= data.max.z; k++) {
00193     for (int j = data.min.y; j <= data.max.y; j++) {
00194       mac->u(data.min.x, j, k) = src->value(data.min.x, j, k).x;
00195       mac->u(data.max.x + 1, j, k) = src->value(data.max.x, j, k).x;
00196     }
00197   }
00198 
00199   // Set the v edge value to their closest voxel
00200   for (int k = data.min.z; k <= data.max.z; k++) {
00201     for (int i = data.min.x; i <= data.max.x; i++) {
00202       mac->v(i, data.min.y, k) = src->value(i, data.min.y, k).y;
00203       mac->v(i, data.max.y + 1, k) = src->value(i, data.max.y, k).y;
00204     }
00205   }
00206 
00207   // Set the w edge value to their closest voxel
00208   for (int j = data.min.y; j <= data.max.y; j++) {
00209     for (int i = data.min.x; i <= data.max.x; i++) {
00210       mac->w(i, j, data.min.z) = src->value(i, j, data.min.z).z;
00211       mac->w(i, j, data.max.z + 1) = src->value(i, j, data.max.z).z;
00212     }
00213   }
00214 
00215   // Loop over internal u values
00216   for (int k = data.min.z; k <= data.max.z; ++k) {
00217     for (int j = data.min.y; j <= data.max.y; ++j) {
00218       for (int i = data.min.x + 1; i <= data.max.x; ++i) {
00219         mac->u(i, j, k) = 
00220           (src->value(i, j, k).x + src->value(i - 1, j, k).x) * 0.5;
00221       }
00222     }
00223   }
00224 
00225   // Loop over internal v values
00226   for (int k = data.min.z; k <= data.max.z; ++k) {
00227     for (int j = data.min.y + 1; j <= data.max.y; ++j) {
00228       for (int i = data.min.x; i <= data.max.x; ++i) {
00229         mac->v(i, j, k) = 
00230           (src->value(i, j, k).y + src->value(i, j - 1, k).y) * 0.5;
00231       }
00232     }
00233   }
00234 
00235   // Loop over internal w values
00236   for (int k = data.min.z + 1; k <= data.max.z; ++k) {
00237     for (int j = data.min.y; j <= data.max.y; ++j) {
00238       for (int i = data.min.x; i <= data.max.x; ++i) {
00239         mac->w(i, j, k) = 
00240           (src->value(i, j, k).z + src->value(i, j, k - 1).z) * 0.5;
00241       }
00242     }
00243   }
00244 }
00245 
00246 //----------------------------------------------------------------------------//
00247 
00248 FIELD3D_NAMESPACE_HEADER_CLOSE
00249 
00250 //----------------------------------------------------------------------------//
00251 
00252 #endif // Include guard