00001 #include "SundanceTetQuadrature.hpp" 00002 #include "PlayaExceptions.hpp" 00003 #include "SundanceOut.hpp" 00004 00005 using namespace Sundance; 00006 using namespace Sundance; 00007 using namespace Sundance; 00008 using namespace Sundance; 00009 using namespace Sundance; 00010 using namespace Teuchos; 00011 00012 00013 void TetQuadrature::getPoints(int order, Array<double>& wgt, 00014 Array<double>& x, 00015 Array<double>& y, 00016 Array<double>& z) 00017 { 00018 Array<double> w; 00019 Array<int> multiplicity; 00020 Array<Array<double> > q; 00021 00022 if (order==1) 00023 { 00024 multiplicity = tuple(1); 00025 w = tuple(1.0); 00026 q.resize(1); 00027 q[0] = tuple(0.25); 00028 } 00029 else if (order==2) 00030 { 00031 multiplicity = tuple(4); 00032 w = tuple(0.25); 00033 q.resize(1); 00034 q[0] = tuple(0.5854101966249685, 0.1381966011250105); 00035 } 00036 else if (order==4) 00037 { 00038 multiplicity = tuple(4, 12); 00039 w = tuple(0.05037379410012282, 0.06654206863329239); 00040 q.resize(2); 00041 q[0] = tuple(0.7716429020672371, 0.7611903264425430e-01); 00042 q[1] = tuple(0.1197005277978019, 0.7183164526766925e-01, 0.4042339134672644); 00043 } 00044 else if (order==6) 00045 { 00046 multiplicity = tuple(1, 4, 12, 12); 00047 w = tuple(0.9040129046014750e-01, 0.1911983427899124e-01, 00048 0.4361493840666568e-01, 0.2581167596199161e-01); 00049 q.resize(4); 00050 q[0] = tuple(0.25); 00051 q[1] = tuple(0.8277192480479295, 0.5742691731735683e-01); 00052 q[2] = tuple(0.5135188412556341e-01, 0.4860510285706072, 0.2312985436519147); 00053 q[3] = tuple(0.2967538129690260, 0.6081079894015281, 0.4756909881472290e-01); 00054 } 00055 else 00056 { 00057 #ifndef TRILINOS_7 00058 SUNDANCE_ERROR("symmetric quadrature rule order " 00059 << order << 00060 " not available for triangles"); 00061 #else 00062 SUNDANCE_ERROR7("symmetric quadrature rule order " 00063 << order << 00064 " not available for triangles"); 00065 #endif 00066 } 00067 00068 for (int i=0; i<q.length(); i++) 00069 { 00070 Array<Array<double> > qPerm; 00071 permute(multiplicity[i], q[i], qPerm); 00072 for (int j=0; j<multiplicity[i]; j++) 00073 { 00074 x.append(qPerm[j][0]); 00075 y.append(qPerm[j][1]); 00076 z.append(qPerm[j][2]); 00077 wgt.append(w[i]); 00078 } 00079 } 00080 } 00081 00082 bool TetQuadrature::supportsOrder(int order) 00083 { 00084 if (order==1 || order==2 || order==4 || order==6) return true; 00085 return false; 00086 } 00087 00088 void TetQuadrature::permute(int m, const Array<double>& q, 00089 Array<Array<double> >& qPerm) 00090 { 00091 qPerm.resize(m); 00092 if (m==1) 00093 { 00094 qPerm[0] = tuple(q[0], q[0], q[0], q[0]); 00095 } 00096 else if (m==4) 00097 { 00098 qPerm[0] = tuple(q[0], q[1], q[1], q[1]); 00099 qPerm[1] = tuple(q[1], q[0], q[1], q[1]); 00100 qPerm[2] = tuple(q[1], q[1], q[0], q[1]); 00101 qPerm[3] = tuple(q[1], q[1], q[1], q[0]); 00102 } 00103 else if (m==12) 00104 { 00105 qPerm[0] = tuple(q[0], q[1], q[2], q[2]); 00106 qPerm[1] = tuple(q[0], q[2], q[1], q[2]); 00107 qPerm[2] = tuple(q[0], q[2], q[2], q[1]); 00108 qPerm[3] = tuple(q[1], q[0], q[2], q[2]); 00109 qPerm[4] = tuple(q[2], q[0], q[1], q[2]); 00110 qPerm[5] = tuple(q[2], q[0], q[2], q[1]); 00111 qPerm[6] = tuple(q[1], q[2], q[0], q[2]); 00112 qPerm[7] = tuple(q[2], q[1], q[0], q[2]); 00113 qPerm[8] = tuple(q[2], q[2], q[0], q[1]); 00114 qPerm[9] = tuple(q[1], q[2], q[2], q[0]); 00115 qPerm[10] = tuple(q[2], q[1], q[2], q[0]); 00116 qPerm[11] = tuple(q[2], q[2], q[1], q[0]); 00117 } 00118 else 00119 { 00120 #ifndef TRILINOS_7 00121 SUNDANCE_ERROR("invalid multiplicity " 00122 << m << 00123 " in TetQuadrature::permute()"); 00124 #else 00125 SUNDANCE_ERROR7("invalid multiplicity " 00126 << m << 00127 " in TetQuadrature::permute()"); 00128 #endif 00129 } 00130 } 00131 00132 bool TetQuadrature::test(int p) 00133 { 00134 Array<double> w; 00135 Array<double> x; 00136 Array<double> y; 00137 Array<double> z; 00138 00139 getPoints(p, w, x, y, z); 00140 bool pass = true; 00141 00142 for (int a=0; a<=p; a++) 00143 { 00144 for (int b=0; b<p-a; b++) 00145 { 00146 int cMax = p - a - b; 00147 for (int c=0; c<=cMax; c++) 00148 { 00149 int dMax = p - a - b - c; 00150 for (int d=0; d<=dMax; d++) 00151 { 00152 double sum = 0.0; 00153 for (int q=0; q<w.length(); q++) 00154 { 00155 sum += (1.0/6.0)*w[q] * pow(x[q], (double) a) * pow(y[q], (double) b) 00156 * pow(z[q], (double) c) 00157 * pow(1.0 - x[q] - y[q] - z[q], (double) d); 00158 } 00159 double err = fabs(sum - exact(a,b,c,d)); 00160 bool localPass = err < 1.0e-14; 00161 pass = pass && localPass; 00162 if (!localPass) 00163 { 00164 fprintf(stderr, "order=%d m (%d, %d, %d %d) q=%22.15g exact=%22.15g\n", p, a, b, c, d, sum, exact(a, b, c, d)); 00165 std::cerr << "error = " << err << std::endl; 00166 } 00167 } 00168 } 00169 } 00170 } 00171 return pass; 00172 } 00173 00174 double TetQuadrature::exact(int a, int b, int c, int d) 00175 { 00176 return fact(a)*fact(b)*fact(c)*fact(d)/fact(a+b+c+d+3); 00177 } 00178 00179 double TetQuadrature::fact(int x) 00180 { 00181 if (x==0) return 1.0; 00182 return x*fact(x-1); 00183 } 00184