00001 #include "SundanceTriangleQuadrature.hpp" 00002 #include "PlayaExceptions.hpp" 00003 #include "SundanceOut.hpp" 00004 #include "SundanceGauss1D.hpp" 00005 #include "PlayaTabs.hpp" 00006 00007 using namespace Sundance; 00008 using namespace Sundance; 00009 using namespace Sundance; 00010 using namespace Sundance; 00011 using namespace Sundance; 00012 using namespace Teuchos; 00013 00014 void TriangleQuadrature::getPoints(int order, Array<double>& wgt, 00015 Array<double>& x, 00016 Array<double>& y) 00017 { 00018 if (!getSymmetricPoints(order, wgt, x, y)) 00019 { 00020 getNonsymmetricPoints(order, wgt, x, y); 00021 } 00022 } 00023 00024 00025 bool TriangleQuadrature::getSymmetricPoints(int order, Array<double>& wgt, 00026 Array<double>& x, 00027 Array<double>& y) 00028 { 00029 00030 int np; 00031 Array<double> w; 00032 Array<int> multiplicity; 00033 Array<Array<double> > q; 00034 00035 if (order==1) 00036 { 00037 multiplicity = tuple(1); 00038 np = 1; 00039 w = tuple(1.0); 00040 q.resize(1); 00041 q[0] = tuple(1.0/3.0, 1.0/3.0, 1.0/3.0); 00042 } 00043 else if (order==2) 00044 { 00045 multiplicity = tuple(3); 00046 np = 3; 00047 w = tuple(1.0/3.0); 00048 q.resize(1); 00049 q[0] = tuple(2.0/3.0, 1.0/6.0, 1.0/6.0); 00050 } 00051 else if (order==3) 00052 { 00053 multiplicity = tuple(6); 00054 np = 6; 00055 w = tuple(1.0/6.0); 00056 q.resize(1); 00057 q[0] = tuple(0.659027622374092, 0.231933368553031, 0.109039009072877); 00058 } 00059 else if (order==4) 00060 { 00061 multiplicity = tuple(3, 3); 00062 np = 6; 00063 w = tuple(0.109951743655322, 0.223381589678011); 00064 q.resize(2); 00065 q[0] = tuple(0.816847572980459, 0.091576213509771, 0.091576213509771); 00066 q[1] = tuple(0.108103018168070, 0.445948490915965, 0.445948490915965); 00067 } 00068 else if (order==5) 00069 { 00070 multiplicity = tuple(1, 3, 3); 00071 np = 7; 00072 q.resize(3); 00073 w = tuple(0.22500000000000, 0.125939180544827, 0.132394152788506); 00074 q[0] = tuple(1.0/3.0, 1.0/3.0, 1.0/3.0); 00075 q[1] = tuple(0.797426985353087, 0.101286507323456, 0.101286507323456); 00076 q[2] = tuple(0.059715871789770, 0.470142064105115, 0.470142064105115); 00077 } 00078 else if (order==6) 00079 { 00080 multiplicity = tuple(3, 3, 6); 00081 np = 12; 00082 q.resize(3); 00083 w = tuple(0.050844906370207, 0.116786275726379, 0.082851075618374); 00084 q[0] = tuple(0.873821971016996, 0.063089014491502, 0.063089014491502); 00085 q[1] = tuple(0.501426509658179, 0.249286745170910, 0.249286745170910); 00086 q[2] = tuple(0.636502499121399, 0.310352451033784, 0.053145049844817); 00087 } 00088 else 00089 { 00090 return false; 00091 } 00092 00093 for (int i=0; i<q.length(); i++) 00094 { 00095 Array<Array<double> > qPerm; 00096 permute(multiplicity[i], q[i], qPerm); 00097 for (int j=0; j<multiplicity[i]; j++) 00098 { 00099 x.append(qPerm[j][0]); 00100 y.append(qPerm[j][1]); 00101 wgt.append(w[i]); 00102 } 00103 } 00104 00105 return true; 00106 } 00107 00108 00109 00110 00111 void TriangleQuadrature::getNonsymmetricPoints(int order, Array<double>& wgt, 00112 Array<double>& x, 00113 Array<double>& y) 00114 { 00115 int nNodes = (order+3)/2; 00116 Gauss1D rule(nNodes, -1.0, 1.0); 00117 Array<double> s = rule.nodes(); 00118 Array<double> t = s; 00119 Array<double> w = rule.weights(); 00120 int n = rule.nPoints(); 00121 00122 wgt.resize(n*n); 00123 x.resize(n*n); 00124 y.resize(n*n); 00125 00126 int k=0; 00127 for (int i=0; i<n; i++) 00128 { 00129 double p = (1.0+s[i])/2.0; 00130 double J = 1.0-p; 00131 for (int j=0; j<n; j++, k++) 00132 { 00133 double q = (1.0 - p)*(1.0+t[j])/2.0; 00134 x[k] = p; 00135 y[k] = q; 00136 wgt[k] = 0.5*w[i]*w[j]*J; 00137 } 00138 } 00139 } 00140 00141 00142 void TriangleQuadrature::permute(int m, const Array<double>& q, 00143 Array<Array<double> >& qPerm) 00144 { 00145 qPerm.resize(m); 00146 if (m==1) 00147 { 00148 qPerm[0] = q; 00149 } 00150 else if (m==3) 00151 { 00152 qPerm[0] = tuple(q[0], q[1], q[2]); 00153 qPerm[1] = tuple(q[1], q[0], q[2]); 00154 qPerm[2] = tuple(q[2], q[1], q[0]); 00155 } 00156 else if (m==6) 00157 { 00158 qPerm[0] = tuple(q[0], q[1], q[2]); 00159 qPerm[1] = tuple(q[0], q[2], q[1]); 00160 qPerm[2] = tuple(q[1], q[0], q[2]); 00161 qPerm[3] = tuple(q[1], q[2], q[0]); 00162 qPerm[4] = tuple(q[2], q[1], q[0]); 00163 qPerm[5] = tuple(q[2], q[0], q[1]); 00164 } 00165 else 00166 { 00167 #ifndef TRILINOS_7 00168 SUNDANCE_ERROR("invalid multiplicity " 00169 << m << 00170 " in TriangleQuadrature::permute()"); 00171 #else 00172 SUNDANCE_ERROR7("invalid multiplicity " 00173 << m << 00174 " in TriangleQuadrature::permute()"); 00175 #endif 00176 } 00177 } 00178 00179 bool TriangleQuadrature::test(int p) 00180 { 00181 Array<double> w; 00182 Array<double> x; 00183 Array<double> y; 00184 00185 getPoints(p, w, x, y); 00186 bool pass = true; 00187 00188 for (int a=0; a<=p; a++) 00189 { 00190 for (int b=0; b<p-a; b++) 00191 { 00192 int cMax = p - a - b; 00193 for (int c=0; c<=cMax; c++) 00194 { 00195 double sum = 0.0; 00196 for (int q=0; q<w.length(); q++) 00197 { 00198 sum += 0.5*w[q] * pow(x[q], (double) a) * pow(y[q], (double) b) 00199 * pow(1.0 - x[q] - y[q], (double) c); 00200 } 00201 double err = fabs(sum - exact(a,b,c)); 00202 bool localPass = err < 1.0e-14; 00203 pass = pass && localPass; 00204 if (!localPass) 00205 { 00206 fprintf(stderr, "order=%d m (%d, %d, %d) q=%22.15g exact=%22.15g\n", p, a, b, c, sum, exact(a, b, c)); 00207 std::cerr << "error = " << err << std::endl; 00208 } 00209 } 00210 } 00211 } 00212 return pass; 00213 } 00214 00215 double TriangleQuadrature::exact(int a, int b, int c) 00216 { 00217 return fact(a)*fact(b)*fact(c)/fact(a+b+c+2); 00218 } 00219 00220 double TriangleQuadrature::fact(int x) 00221 { 00222 if (x==0) return 1.0; 00223 return x*fact(x-1); 00224 } 00225 00226