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 "SundancePointwiseUserDefFunctor.hpp"
00032 #include "PlayaTabs.hpp"
00033 #include "SundanceOut.hpp"
00034
00035 using namespace Sundance;
00036 using namespace Sundance;
00037
00038 using namespace Sundance;
00039 using namespace Teuchos;
00040
00041
00042 PointwiseUserDefFunctor0::PointwiseUserDefFunctor0(const std::string& name,
00043 int domainDim,
00044 int rangeDim)
00045 : UserDefFunctor(name, domainDim, rangeDim)
00046 {}
00047
00048 void PointwiseUserDefFunctor0::evaluationCallback(int nPoints, int maxDiffOrder,
00049 const double** in,
00050 double** out) const
00051 {
00052 TEUCHOS_TEST_FOR_EXCEPTION(maxDiffOrder > 0, std::runtime_error,
00053 "diff order = " << maxDiffOrder
00054 << " not supported for functor "
00055 << name());
00056
00057 static Array<double> x;
00058 x.resize(domainDim());
00059
00060 static Array<double> f;
00061 f.resize(rangeDim());
00062
00063 for (int i=0; i<nPoints; i++)
00064 {
00065 for (int j=0; j<domainDim(); j++) x[j] = in[j][i];
00066 const double* xp = &(x[0]);
00067 double* fp = &(f[0]);
00068 eval0(xp, fp);
00069 for (int j=0; j<rangeDim(); j++) out[j][i] = fp[j];
00070 }
00071 }
00072
00073 PointwiseUserDefFunctor1::PointwiseUserDefFunctor1(const std::string& name,
00074 int domainDim,
00075 int rangeDim)
00076 : PointwiseUserDefFunctor0(name, domainDim, rangeDim)
00077 {}
00078
00079 void PointwiseUserDefFunctor1::evaluationCallback(int nPoints, int maxDiffOrder,
00080 const double** in,
00081 double** out) const
00082 {
00083 TEUCHOS_TEST_FOR_EXCEPTION(maxDiffOrder > 1, std::runtime_error,
00084 "diff order = " << maxDiffOrder
00085 << " not supported for functor "
00086 << name());
00087
00088 static Array<double> x;
00089 x.resize(domainDim());
00090
00091 static Array<double> f;
00092 if (maxDiffOrder==1) f.resize(rangeDim() * (1 + domainDim()) );
00093 else f.resize(rangeDim());
00094 double* fp = &(f[0]);
00095
00096 for (int i=0; i<nPoints; i++)
00097 {
00098 for (int j=0; j<domainDim(); j++) x[j] = in[j][i];
00099 const double* xp = &(x[0]);
00100
00101 if (maxDiffOrder==1)
00102 {
00103 double* dfp = &(f[rangeDim()]);
00104 eval1(xp, fp, dfp);
00105 }
00106 else eval0(xp, fp);
00107 for (int j=0; j<f.size(); j++) out[j][i] = fp[j];
00108 }
00109 }
00110
00111
00112 void PointwiseUserDefFunctor1::eval0(const double* in, double* out) const
00113 {
00114 static Array<double> dummy;
00115 dummy.resize(domainDim() * rangeDim());
00116
00117 eval1(in, out, &(dummy[0]));
00118 }
00119
00120
00121
00122 PointwiseUserDefFunctor2::PointwiseUserDefFunctor2(const std::string& name,
00123 int domainDim,
00124 int rangeDim)
00125 : PointwiseUserDefFunctor1(name, domainDim, rangeDim)
00126 {}
00127
00128 void PointwiseUserDefFunctor2::evaluationCallback(int nPoints, int maxDiffOrder,
00129 const double** in,
00130 double** out) const
00131 {
00132 TEUCHOS_TEST_FOR_EXCEPTION(maxDiffOrder > 2 || maxDiffOrder < 0, std::runtime_error,
00133 "diff order = " << maxDiffOrder
00134 << " not supported for functor "
00135 << name());
00136
00137 int nTotal = 1;
00138 int numFirst = domainDim();
00139 int numSecond = domainDim()*(domainDim()+1)/2;
00140
00141 static Array<double> x;
00142 x.resize(domainDim());
00143
00144 static Array<double> f;
00145
00146 if (maxDiffOrder > 0) nTotal += numFirst;
00147 if (maxDiffOrder > 1) nTotal += numSecond;
00148
00149 f.resize(rangeDim() * nTotal);
00150
00151 double* fp = &(f[0]);
00152
00153 for (int i=0; i<nPoints; i++)
00154 {
00155 for (int j=0; j<domainDim(); j++) x[j] = in[j][i];
00156 const double* xp = &(x[0]);
00157
00158 if (maxDiffOrder==0)
00159 {
00160 eval0(xp, fp);
00161 }
00162 else if (maxDiffOrder==1)
00163 {
00164 double* dfp = &(f[rangeDim()]);
00165 eval1(xp, fp, dfp);
00166 }
00167 else if (maxDiffOrder==2)
00168 {
00169 double* dfp = &(f[rangeDim()]);
00170 double* d2fp = &(f[rangeDim()*(1 + domainDim())]);
00171 eval2(xp, fp, dfp, d2fp);
00172 }
00173 else
00174 {
00175 TEUCHOS_TEST_FOR_EXCEPT(true);
00176 }
00177 for (int j=0; j<f.size(); j++) out[j][i] = fp[j];
00178 }
00179 }
00180
00181
00182 void PointwiseUserDefFunctor2::eval0(const double* in, double* f) const
00183 {
00184 static Array<double> dummy1;
00185 static Array<double> dummy2;
00186 dummy1.resize(rangeDim() * domainDim());
00187 dummy2.resize(rangeDim() * domainDim()*(domainDim()+1)/2);
00188
00189 eval2(in, f, &(dummy1[0]), &(dummy2[0]));
00190 }
00191
00192
00193 void PointwiseUserDefFunctor2::eval1(const double* in, double* f, double* df) const
00194 {
00195 static Array<double> dummy2;
00196 dummy2.resize(rangeDim() * domainDim()*(domainDim()+1)/2);
00197
00198 eval2(in, f, df, &(dummy2[0]));
00199 }