00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceFourier.hpp"
00032 #include "SundanceADReal.hpp"
00033 #include "PlayaExceptions.hpp"
00034 #include "SundanceSpatialDerivSpecifier.hpp"
00035 #include "SundancePoint.hpp"
00036 #include "SundanceObjectWithVerbosity.hpp"
00037 #include "SundanceOut.hpp"
00038
00039 using namespace Sundance;
00040 using namespace Teuchos;
00041 using std::endl;
00042
00043
00044 Fourier::Fourier(int N) : N_(N)
00045 {}
00046
00047 bool Fourier::supportsCellTypePair(
00048 const CellType& maximalCellType,
00049 const CellType& cellType
00050 ) const
00051 {
00052 switch(maximalCellType)
00053 {
00054 case LineCell:
00055 switch(cellType)
00056 {
00057 case LineCell:
00058 case PointCell:
00059 return true;
00060 default:
00061 return false;
00062 }
00063 default:
00064 return false;
00065 }
00066 }
00067
00068 void Fourier::print(std::ostream& os) const
00069 {
00070 os << "Fourier(" << N_ << ")";
00071 }
00072
00073 int Fourier::nReferenceDOFsWithoutFacets(
00074 const CellType& maximalCellType,
00075 const CellType& cellType
00076 ) const
00077 {
00078 switch(cellType)
00079 {
00080 case PointCell:
00081 return 0;
00082 case LineCell:
00083 return 2*N_+1;
00084 default:
00085 TEUCHOS_TEST_FOR_EXCEPTION( true , std::invalid_argument , "illegal combination of cell type and maximal cell type" );
00086 return -1;
00087 }
00088 }
00089
00090 Array<int> Fourier::makeRange(int low, int high)
00091 {
00092 if (high < low) return Array<int>();
00093
00094 Array<int> rtn(high-low+1);
00095 for (int i=0; i<rtn.length(); i++) rtn[i] = low+i;
00096 return rtn;
00097 }
00098
00099 void Fourier::getReferenceDOFs(
00100 const CellType& maximalCellType,
00101 const CellType& cellType,
00102 Array<Array<Array<int> > >& dofs) const
00103 {
00104 typedef Array<int> Aint;
00105
00106 switch(cellType)
00107 {
00108 case LineCell:
00109 {
00110 dofs.resize(2);
00111 dofs[0] = Array<Array<int> >();
00112 dofs[1] = tuple(makeRange(0, 2*N_));
00113 }
00114 break;
00115 default:
00116 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00117 << cellType << " not implemented in Fourier basis");
00118 }
00119 }
00120
00121
00122 void Fourier::refEval(
00123 const CellType& cellType,
00124 const Array<Point>& pts,
00125 const SpatialDerivSpecifier& sds,
00126 Array<Array<Array<double> > >& result,
00127 int verbosity) const
00128 {
00129 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00130 std::runtime_error,
00131 "cannot evaluate spatial derivative " << sds << " on Fourier basis");
00132 const MultiIndex& deriv = sds.mi();
00133 typedef Array<double> Adouble;
00134 result.resize(1);
00135 result[0].resize(pts.length());
00136
00137 switch(cellType)
00138 {
00139 case LineCell:
00140 for (int i=0; i<pts.length(); i++)
00141 {
00142 evalOnLine(pts[i], deriv, result[0][i]);
00143 }
00144 break;
00145 default:
00146 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00147 "Fourier::refEval() unimplemented for cell type "
00148 << cellType);
00149 }
00150 }
00151
00152
00153
00154 void Fourier::evalOnLine(const Point& pt,
00155 const MultiIndex& deriv,
00156 Array<double>& result) const
00157 {
00158 result.resize(2*N_+1);
00159 double x = pt[0];
00160 const double pi = 4.0*atan(1.0);
00161
00162 if (deriv.order()==0)
00163 {
00164 result[0] = 1.0;
00165 for (int n=1; n<=N_; n++)
00166 {
00167 result[2*n-1]=::cos(2.0*pi*n*x);
00168 result[2*n]=::sin(2.0*pi*n*x);
00169 }
00170 }
00171 else if (deriv.order()==1)
00172 {
00173 result[0] = 0.0;
00174 for (int n=1; n<=N_; n++)
00175 {
00176 result[2*n-1]=-2.0*pi*n*::sin(2.0*pi*n*x);
00177 result[2*n]=2.0*n*pi*::cos(2.0*pi*n*x);
00178 }
00179 }
00180 else
00181 {
00182 TEUCHOS_TEST_FOR_EXCEPT(true);
00183 }
00184
00185 Out::os() << "quad point=" << x << endl;
00186 Out::os() << "bas: " << result << endl;
00187 }
00188