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 "SundanceRaviartThomas.hpp"
00032 #include "SundancePoint.hpp"
00033 #include "SundanceCellType.hpp"
00034 #include "SundanceSpatialDerivSpecifier.hpp"
00035 #include "PlayaExceptions.hpp"
00036 #include "SundanceTypeUtils.hpp"
00037 #include "SundanceObjectWithVerbosity.hpp"
00038 #include "SundanceOut.hpp"
00039 #include "SundanceADReal.hpp"
00040
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043
00044
00045 RaviartThomas::RaviartThomas(int spatialDim)
00046 : HDivVectorBasis(spatialDim)
00047 {}
00048
00049 bool RaviartThomas::supportsCellTypePair(
00050 const CellType& maximalCellType,
00051 const CellType& cellType
00052 ) const
00053 {
00054 switch(maximalCellType)
00055 {
00056 case TriangleCell:
00057 switch(cellType)
00058 {
00059 case TriangleCell:
00060 case LineCell:
00061 case PointCell:
00062 return true;
00063 default:
00064 return false;
00065 }
00066 case TetCell:
00067 switch(cellType)
00068 {
00069 case TetCell:
00070 case TriangleCell:
00071 case LineCell:
00072 case PointCell:
00073 return true;
00074 default:
00075 return false;
00076 }
00077 default:
00078 return false;
00079 }
00080 }
00081
00082 std::string RaviartThomas::description() const
00083 {
00084 return "RaviartThomas()";
00085 }
00086
00087 int RaviartThomas::nReferenceDOFsWithoutFacets(
00088 const CellType& maximalCellType,
00089 const CellType& cellType
00090 ) const
00091 {
00092 if (dimension(cellType) == dimension(maximalCellType)-1)
00093 {
00094 return 1;
00095 }
00096 else
00097 {
00098 return 0;
00099 }
00100 }
00101
00102 void RaviartThomas::getReferenceDOFs(
00103 const CellType& maximalCellType,
00104 const CellType& cellType,
00105 Array<Array<Array<int> > >& dofs) const
00106 {
00107 switch(cellType)
00108 {
00109 case PointCell:
00110 dofs.resize(1);
00111 dofs[0] = tuple(Array<int>());
00112 return;
00113 case LineCell:
00114 dofs.resize(2);
00115 dofs[0] = tuple(Array<int>());
00116 dofs[1] = tuple<Array<int> >(tuple(0));
00117 return;
00118 case TriangleCell:
00119 dofs.resize(3);
00120 dofs[0] = tuple(Array<int>());
00121 dofs[1] = tuple<Array<int> >(tuple(0), tuple(1), tuple(2));
00122 dofs[2] = tuple(Array<int>());
00123 return;
00124 case TetCell:
00125 dofs.resize(4);
00126 dofs[0] = tuple(Array<int>());
00127 dofs[1] = tuple<Array<int> >(tuple(0), tuple(1), tuple(2), tuple(3));
00128 dofs[2] = tuple(Array<int>());
00129 dofs[3] = tuple(Array<int>());
00130 return;
00131 default:
00132 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Cell type "
00133 << cellType << " not implemented in RaviartThomas basis");
00134 }
00135 }
00136
00137
00138
00139 bool RaviartThomas::lessThan(const BasisDOFTopologyBase* other) const
00140 {
00141 if (typeLessThan(this, other)) return true;
00142 if (typeLessThan(other, this)) return false;
00143
00144 return false;
00145 }
00146
00147
00148 void RaviartThomas::refEval(
00149 const CellType& cellType,
00150 const Array<Point>& pts,
00151 const SpatialDerivSpecifier& sds,
00152 Array<Array<Array<double> > >& result,
00153 int verbosity) const
00154 {
00155 const MultiIndex& deriv = sds.mi();
00156 switch(cellType) {
00157 case PointCell:
00158 {
00159 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "evaluation of RaviartThomas elements for PointCell not supported");
00160 }
00161 break;
00162 case LineCell:
00163 {
00164 result.resize(1);
00165 result[0].resize(pts.length());
00166 Array<ADReal> tmp;
00167 tmp.resize(2);
00168 for (int i=0;i<pts.length();i++) {
00169 ADReal x = ADReal( pts[i][0] , 0 , 1 );
00170 ADReal one(1.0,1);
00171 result[0][i].resize(2);
00172 tmp[0] = one - x;
00173 tmp[1] = x;
00174 if (deriv.order()==0) {
00175 for (int j=0;j<tmp.length();j++) {
00176 result[0][i][j] = tmp[j].value();
00177 }
00178 }
00179 else {
00180 for (int j=0;j<tmp.length();j++) {
00181 result[0][i][j] = tmp[j].gradient()[deriv.firstOrderDirection()];
00182 }
00183 }
00184 }
00185 }
00186 break;
00187 case TriangleCell:
00188 {
00189 result.resize(2);
00190 result[0].resize(pts.length());
00191 result[1].resize(pts.length());
00192 Array<ADReal> tmp0;
00193 Array<ADReal> tmp1;
00194 tmp0.resize(3);
00195 tmp1.resize(3);
00196 for (int i=0;i<pts.length();i++) {
00197 ADReal x = ADReal(pts[i][0],0,2);
00198 ADReal y = ADReal(pts[i][1],1,2);
00199 ADReal one(1.0,2);
00200 ADReal rt2(std::sqrt(2.0),2);
00201 result[0][i].resize(3);
00202 result[1][i].resize(3);
00203 tmp0[0] = x;
00204 tmp1[0] = y - one;
00205 tmp0[1] = rt2 * x;
00206 tmp1[1] = rt2 * y;
00207 tmp0[2] = x - one;
00208 tmp1[2] = y;
00209 if (deriv.order()==0) {
00210 for (int j=0;j<tmp0.length();j++) {
00211 result[0][i][j] = tmp0[j].value();
00212 result[1][i][j] = tmp1[j].value();
00213 }
00214 }
00215 else {
00216 for (int j=0;j<tmp0.length();j++) {
00217 result[0][i][j] = tmp0[j].gradient()[deriv.firstOrderDirection()];
00218 result[1][i][j] = tmp1[j].gradient()[deriv.firstOrderDirection()];
00219 }
00220 }
00221 }
00222 }
00223 break;
00224 case TetCell:
00225 {
00226 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "evaluation of RaviartThomas elements for TetCell not supported");
00227 }
00228 break;
00229 default:
00230 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "evaluation of RaviartThomas elements unknown cell type");
00231 break;
00232 }
00233
00234 return;
00235 }
00236
00237
00238 void RaviartThomas::print(std::ostream& os) const
00239 {
00240 os << "RaviartThomas(" << dim() << ")";
00241 }