PartialDerivative.hpp

00001 #ifndef __FIBER_VECTOR_PARTIAL_DERIVATIVE_HPP
00002 #define __FIBER_VECTOR_PARTIAL_DERIVATIVE_HPP
00003 
00004 #include <eagle/Vector.hpp>
00005 #include "MultiArray.hpp"
00006 
00007 namespace Fiber
00008 {
00009 
00010 
00016 template <int IndexChart, typename Value, typename ResultType>
00017 Eagle::Vector<ResultType, IndexChart>&
00018 PartialDerivative(Eagle::Vector<ResultType, IndexChart>&result, const MultiArray<IndexChart, Value>&ValueField, const MultiIndex<IndexChart>&I)
00019 {
00020 
00021 const MultiIndex<IndexChart>&Sz = ValueField.Size();
00022 
00023         for(int i=0; i<IndexChart; i++)
00024         {
00025         MultiIndex<IndexChart> DiffI = MultiIndex<IndexChart>::Axis(i); 
00026 
00027                 if (I[i]<1) // on left boundary, need to do forward-differencing
00028                 {
00029                         result[i] = ValueField[ I + DiffI ] - ValueField[ I ]; 
00030                 } 
00031                 else if (I[i]>=Sz[i]-1) // on right boundary, need to do backward-differencing
00032                 {
00033                         result[i] = ValueField[ I ] - ValueField[ I - DiffI ]; 
00034                 }
00035                 else
00036                 {
00037                         result[i] =  0.5*(ValueField[ I + DiffI ] - ValueField[ I - DiffI ]);
00038                 }
00039         }
00040 
00041         return result;
00042 }
00043 
00044 
00045 } /* namespace Fiber */ 
00046 
00047 #endif /* __FIBER_VECTOR_PARTIAL_DERIVATIVE_HPP */
00048