00001 #ifndef OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_DIV_GRAD_H
00002 #define OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_DIV_GRAD_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 #include <OpenTissue/core/containers/grid/util/grid_gradient.h>
00013 #include <OpenTissue/core/containers/grid/util/grid_coord2idx.h>
00014 #include <cmath>
00015
00016 namespace OpenTissue
00017 {
00018 namespace grid
00019 {
00020
00028 template<typename grid_type>
00029 inline void div_grad(
00030 grid_type const & phi
00031 , grid_type & M
00032 )
00033 {
00034 using std::min;
00035
00036 typedef OpenTissue::math::Vector3< typename grid_type::value_type> vector3_type;
00037
00038 typedef typename grid_type::iterator iterator;
00039 typedef typename grid_type::const_index_iterator const_index_iterator;
00040 typedef typename grid_type::value_type value_type;
00041 typedef typename grid_type::math_types math_types;
00042 typedef typename math_types::real_type real_type;
00043
00044 size_t const & I = phi.I();
00045 size_t const & J = phi.J();
00046 size_t const & K = phi.K();
00047 M.create(phi.min_coord(),phi.max_coord(),I,J,K);
00048 static value_type unused = phi.unused();
00049
00050 vector3_type gp00,gm00,g0p0,g0m0,g00p,g00m;
00051 vector3_type gppp,gppm,gpmp,gpmm,gmpp,gmpm,gmmp,gmmm;
00052 vector3_type gpp0, gpm0, gmp0, gmm0, gp0p, gp0m, gm0p, gm0m, g0pp, g0pm, g0mp, g0mm;
00053
00054 static vector3_type np00 = unit(vector3_type( 1, 0, 0));
00055 static vector3_type nm00 = unit(vector3_type(-1, 0, 0));
00056 static vector3_type n0p0 = unit(vector3_type( 0, 1, 0));
00057 static vector3_type n0m0 = unit(vector3_type( 0,-1, 0));
00058 static vector3_type n00p = unit(vector3_type( 0, 0, 1));
00059 static vector3_type n00m = unit(vector3_type( 0, 0,-1));
00060 static vector3_type nppp = unit(vector3_type( 1, 1, 1));
00061 static vector3_type nppm = unit(vector3_type( 1, 1,-1));
00062 static vector3_type npmp = unit(vector3_type( 1,-1, 1));
00063 static vector3_type npmm = unit(vector3_type( 1,-1,-1));
00064 static vector3_type nmpp = unit(vector3_type(-1, 1, 1));
00065 static vector3_type nmpm = unit(vector3_type(-1, 1,-1));
00066 static vector3_type nmmp = unit(vector3_type(-1,-1, 1));
00067 static vector3_type nmmm = unit(vector3_type(-1,-1,-1));
00068 static vector3_type npp0 = unit(vector3_type( 1, 1, 0));
00069 static vector3_type npm0 = unit(vector3_type( 1,-1, 0));
00070 static vector3_type nmp0 = unit(vector3_type(-1, 1, 0));
00071 static vector3_type nmm0 = unit(vector3_type(-1,-1, 0));
00072 static vector3_type np0p = unit(vector3_type( 1, 0, 1));
00073 static vector3_type np0m = unit(vector3_type( 1, 0,-1));
00074 static vector3_type nm0p = unit(vector3_type(-1, 0, 1));
00075 static vector3_type nm0m = unit(vector3_type(-1, 0,-1));
00076 static vector3_type n0pp = unit(vector3_type( 0, 1, 1));
00077 static vector3_type n0pm = unit(vector3_type( 0, 1,-1));
00078 static vector3_type n0mp = unit(vector3_type( 0,-1, 1));
00079 static vector3_type n0mm = unit(vector3_type( 0,-1,-1));
00080
00081
00082 iterator m = M.begin();
00083 const_index_iterator pbegin = phi.begin();
00084 const_index_iterator pend = phi.end();
00085 const_index_iterator p = pbegin;
00086 for(;p!=pend;++p,++m)
00087 {
00088 if((*p)==unused)
00089 {
00090 *m = value_type();
00091 continue;
00092 }
00093
00094 size_t i = p.i();
00095 size_t j = p.j();
00096 size_t k = p.k();
00097 size_t im1 = ( i ) ? i - 1 : 0;
00098 size_t jm1 = ( j ) ? j - 1 : 0;
00099 size_t km1 = ( k ) ? k - 1 : 0;
00100 size_t ip1 = min( i + 1u, I - 1u );
00101 size_t jp1 = min( j + 1u, J - 1u );
00102 size_t kp1 = min( k + 1u, K - 1u );
00103
00104 gradient(phi, ip1, j, k, gp00);
00105 gradient(phi, im1, j, k, gm00);
00106 gradient(phi, i, jp1, k, g0p0);
00107 gradient(phi, i, jm1, k, g0m0);
00108 gradient(phi, i, j, kp1, g00p);
00109 gradient(phi, i, j, km1, g00m);
00110 gradient(phi, ip1, jp1, k, gpp0);
00111 gradient(phi, ip1, jm1, k, gpm0);
00112 gradient(phi, im1, jp1, k, gmp0);
00113 gradient(phi, im1, jm1, k, gmm0);
00114 gradient(phi, ip1, j, kp1, gp0p);
00115 gradient(phi, ip1, j, km1, gp0m);
00116 gradient(phi, im1, j, kp1, gm0p);
00117 gradient(phi, im1, j, km1, gm0m);
00118 gradient(phi, i, jp1, kp1, g0pp);
00119 gradient(phi, i, jp1, km1, g0pm);
00120 gradient(phi, i, jm1, kp1, g0mp);
00121 gradient(phi, i, jm1, km1, g0mm);
00122 gradient(phi, ip1, jp1, kp1, gppp);
00123 gradient(phi, ip1, jp1, km1, gppm);
00124 gradient(phi, ip1, jm1, kp1, gpmp);
00125 gradient(phi, ip1, jm1, km1, gpmm);
00126 gradient(phi, im1, jp1, kp1, gmpp);
00127 gradient(phi, im1, jp1, km1, gmpm);
00128 gradient(phi, ip1, jm1, kp1, gmmp);
00129 gradient(phi, im1, jm1, km1, gmmm);
00130
00131 real_type flux = real_type();
00132 flux += gp00 * np00;
00133 flux += gm00 * nm00;
00134 flux += g0p0 * n0p0;
00135 flux += g0m0 * n0m0;
00136 flux += g00p * n00p;
00137 flux += g00m * n00m;
00138 flux += gpp0 * npp0;
00139 flux += gpm0 * npm0;
00140 flux += gmp0 * nmp0;
00141 flux += gmm0 * nmm0;
00142 flux += gp0p * np0p;
00143 flux += gp0m * np0m;
00144 flux += gm0p * nm0p;
00145 flux += gm0m * nm0m;
00146 flux += g0pp * n0pp;
00147 flux += g0pm * n0pm;
00148 flux += g0mp * n0mp;
00149 flux += g0mm * n0mm;
00150 flux += gppp * nppp;
00151 flux += gppm * nppm;
00152 flux += gpmp * npmp;
00153 flux += gpmm * npmm;
00154 flux += gmpp * nmpp;
00155 flux += gmpm * nmpm;
00156 flux += gmmp * nmmp;
00157 flux += gmmm * nmmm;
00158
00159 (*m) = static_cast<value_type>(flux);
00160 }
00161 }
00162
00163 }
00164 }
00165
00166
00167 #endif
