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 "SundanceEdgeLocalizedBasis.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 Sundance;
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043
00044
00045 EdgeLocalizedBasis::EdgeLocalizedBasis()
00046 {}
00047
00048 bool EdgeLocalizedBasis::supportsCellTypePair(
00049 const CellType& maximalCellType,
00050 const CellType& cellType
00051 ) const
00052 {
00053 switch(maximalCellType)
00054 {
00055 case TriangleCell:
00056 case TetCell:
00057 case QuadCell:
00058 case BrickCell:
00059 switch(cellType)
00060 {
00061 case LineCell:
00062 return true;
00063 default:
00064 return false;
00065 }
00066 default:
00067 return false;
00068 }
00069 }
00070
00071 void EdgeLocalizedBasis::print(std::ostream& os) const
00072 {
00073 os << "EdgeLocalizedBasis()";
00074 }
00075
00076 int EdgeLocalizedBasis::nReferenceDOFsWithoutFacets(
00077 const CellType& maximalCellType,
00078 const CellType& cellType
00079 ) const
00080 {
00081 switch(cellType)
00082 {
00083 case LineCell:
00084 return 1;
00085 default:
00086 return 0;
00087 }
00088 }
00089
00090 void EdgeLocalizedBasis::getReferenceDOFs(
00091 const CellType& maximalCellType,
00092 const CellType& cellType,
00093 Array<Array<Array<int> > >& dofs) const
00094 {
00095 typedef Array<int> Aint;
00096 switch(cellType)
00097 {
00098 case LineCell:
00099 dofs.resize(2);
00100 dofs[0] = Array<Array<int> >();
00101 dofs[1] = tuple<Aint>(tuple<int>(0));
00102 return;
00103 case TriangleCell:
00104 dofs.resize(3);
00105 dofs[0] = tuple(Array<int>());
00106 dofs[1] = tuple<Aint>(tuple(0), tuple(1), tuple(2));
00107 dofs[2] = tuple(Array<int>());
00108 return;
00109 default:
00110 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00111 << cellType << " not implemented in EdgeLocalizedBasis basis");
00112 }
00113 }
00114
00115
00116
00117 void EdgeLocalizedBasis::refEval(
00118 const CellType& cellType,
00119 const Array<Point>& pts,
00120 const SpatialDerivSpecifier& sds,
00121 Array<Array<Array<double> > >& result,
00122 int verbosity) const
00123 {
00124 TEUCHOS_TEST_FOR_EXCEPTION(!(sds.isPartial() || sds.isIdentity()),
00125 std::runtime_error,
00126 "cannot evaluate spatial derivative " << sds << " on EdgeLocalizedBasis basis");
00127 const MultiIndex& deriv = sds.mi();
00128 typedef Array<double> Adouble;
00129 result.resize(1);
00130 result[0].resize(pts.length());
00131
00132 int dim = dimension(cellType);
00133
00134 if (dim==0)
00135 {
00136 result[0] = tuple<Adouble>(tuple(1.0));
00137 }
00138 else if (dim==1)
00139 {
00140 for (int i=0; i<pts.length(); i++)
00141 {
00142 evalOnLine(pts[i], deriv, result[0][i]);
00143 }
00144 }
00145 else if (dim==2)
00146 {
00147 for (int i=0; i<pts.length(); i++)
00148 {
00149 evalOnTriangle(pts[i], deriv, result[0][i]);
00150 }
00151 }
00152 else if (dim==3)
00153 {
00154 for (int i=0; i<pts.length(); i++)
00155 {
00156 evalOnTet(pts[i], deriv, result[0][i]);
00157 }
00158 }
00159 }
00160
00161
00162
00163
00164 void EdgeLocalizedBasis::evalOnLine(const Point& pt,
00165 const MultiIndex& deriv,
00166 Array<double>& result) const
00167 {
00168 ADReal one(1.0, 1);
00169 result.resize(1);
00170 Array<ADReal> tmp(result.length());
00171
00172 tmp[0] = one;
00173
00174 for (int i=0; i<tmp.length(); i++)
00175 {
00176 if (deriv.order()==0) result[i] = tmp[i].value();
00177 else result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00178 }
00179 }
00180
00181 void EdgeLocalizedBasis::evalOnTriangle(const Point& pt,
00182 const MultiIndex& deriv,
00183 Array<double>& result) const
00184 {
00185 ADReal x = ADReal(pt[0], 0, 2);
00186 ADReal y = ADReal(pt[1], 1, 2);
00187 ADReal one(1.0, 2);
00188 ADReal zero(0.0, 2);
00189
00190 Array<ADReal> tmp;
00191
00192 SUNDANCE_OUT(this->verb() > 3, "x=" << x.value() << " y="
00193 << y.value());
00194
00195 result.resize(3);
00196 tmp.resize(3);
00197
00198 bool onEdge2 = std::fabs(pt[1]) < 1.0e-14;
00199 bool onEdge0 = std::fabs(1.0-pt[0]-pt[1]) < 1.0e-14;
00200 bool onEdge1 = std::fabs(pt[0]) < 1.0e-14;
00201
00202 TEUCHOS_TEST_FOR_EXCEPTION(!(onEdge0 || onEdge1 || onEdge2),
00203 std::runtime_error,
00204 "EdgeLocalizedBasis should not be evaluated at points not on edges");
00205
00206 TEUCHOS_TEST_FOR_EXCEPTION((onEdge0 && onEdge1) || (onEdge1 && onEdge2)
00207 || (onEdge2 && onEdge0), std::runtime_error,
00208 "Ambiguous edge in EdgeLocalizedBasis::evalOnTriangle()");
00209
00210 if (onEdge0)
00211 {
00212 tmp[0] = one;
00213 tmp[1] = zero;
00214 tmp[2] = zero;
00215 }
00216 if (onEdge1)
00217 {
00218 tmp[0] = zero;
00219 tmp[1] = one;
00220 tmp[2] = zero;
00221 }
00222 if (onEdge2)
00223 {
00224 tmp[0] = zero;
00225 tmp[1] = zero;
00226 tmp[2] = one;
00227 }
00228
00229
00230 for (int i=0; i<tmp.length(); i++)
00231 {
00232 SUNDANCE_OUT(this->verb() > 3,
00233 "tmp[" << i << "]=" << tmp[i].value()
00234 << " grad=" << tmp[i].gradient());
00235 if (deriv.order()==0) result[i] = tmp[i].value();
00236 else
00237 result[i] = tmp[i].gradient()[deriv.firstOrderDirection()];
00238 }
00239
00240 }
00241
00242
00243
00244
00245 void EdgeLocalizedBasis::evalOnTet(const Point& pt,
00246 const MultiIndex& deriv,
00247 Array<double>& result) const
00248 {
00249 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00250 "EdgeLocalizedBasis::evalOnTet not implemented");
00251 }