interpolations3D.hh

00001 /* This file is part of the Palabos library.
00002  *
00003  * Copyright (C) 2011 FlowKit Sarl
00004  * Avenue de Chailly 23
00005  * 1012 Lausanne, Switzerland
00006  * E-mail contact: contact@flowkit.com
00007  *
00008  * The most recent release of Palabos can be downloaded at 
00009  * <http://www.palabos.org/>
00010  *
00011  * The library Palabos is free software: you can redistribute it and/or
00012  * modify it under the terms of the GNU Affero General Public License as
00013  * published by the Free Software Foundation, either version 3 of the
00014  * License, or (at your option) any later version.
00015  *
00016  * The library is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU Affero General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU Affero General Public License
00022  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00023 */
00024 
00025 #ifndef INTERPOLATIONS_3D_HH
00026 #define INTERPOLATIONS_3D_HH
00027 
00028 #include "core/globalDefs.h"
00029 #include "core/util.h"
00030 #include "finiteDifference/interpolations3D.h"
00031 #include <vector>
00032 
00033 namespace plb {
00034 
00035 /* ******** Function linearInterpolationCoefficients ********************* */
00036 
00037 template<typename T>
00038 void linearInterpolationCoefficients (
00039         AtomicBlock3D const& block, Array<T,3> const& position,
00040         std::vector<Dot3D>& cellPos, std::vector<T>& weights )
00041 {
00042     cellPos.resize(8);
00043     cellPos[0] = Dot3D( (plint) position[0],         (plint) position[1],         (plint) position[2] );
00044     cellPos[1] = Dot3D( (plint) position[0],         (plint) position[1],         (plint)(position[2]+(T)1.0) );
00045     cellPos[2] = Dot3D( (plint) position[0],         (plint)(position[1]+(T)1.0), (plint) position[2] );
00046     cellPos[3] = Dot3D( (plint) position[0],         (plint)(position[1]+(T)1.0), (plint)(position[2]+(T)1.0) );
00047     cellPos[4] = Dot3D( (plint)(position[0]+(T)1.0), (plint) position[1],         (plint) position[2] );
00048     cellPos[5] = Dot3D( (plint)(position[0]+(T)1.0), (plint) position[1],         (plint)(position[2]+(T)1.0) );
00049     cellPos[6] = Dot3D( (plint)(position[0]+(T)1.0), (plint)(position[1]+(T)1.0), (plint) position[2] );
00050     cellPos[7] = Dot3D( (plint)(position[0]+(T)1.0), (plint)(position[1]+(T)1.0), (plint)(position[2]+(T)1.0) );
00051 
00052     T u = position[0] - (T)cellPos[0].x;
00053     T v = position[1] - (T)cellPos[0].y;
00054     T w = position[2] - (T)cellPos[0].z;
00055 
00056     weights.resize(8);
00057     weights[0] = (1.-u) * (1.-v) * (1.-w);
00058     weights[1] = (1.-u) * (1.-v) * (   w);
00059     weights[2] = (1.-u) * (   v) * (1.-w);
00060     weights[3] = (1.-u) * (   v) * (   w);
00061     weights[4] = (   u) * (1.-v) * (1.-w);
00062     weights[5] = (   u) * (1.-v) * (   w);
00063     weights[6] = (   u) * (   v) * (1.-w);
00064     weights[7] = (   u) * (   v) * (   w);
00065 
00066     PLB_ASSERT( util::fpequal(
00067                         weights[0] + weights[1] + weights[2] +
00068                         weights[3] + weights[4] + weights[5] +
00069                         weights[6] + weights[7],
00070                     (T)1., (T)1.e-6) );
00071 
00072     // Convert cell position to local coordinates.
00073     for (plint iPos=0; iPos<8; ++iPos) {
00074         cellPos[iPos] -= block.getLocation();
00075     }
00076 }
00077 
00078 }  // namespace plb
00079 
00080 #endif  // INTERPOLATIONS_3D_HH