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 "SundanceGaussianQuadrature.hpp"
00032 #include "SundanceGauss1D.hpp"
00033 #include "SundanceTriangleQuadrature.hpp"
00034 #include "SundanceQuadQuadrature.hpp"
00035 #include "SundanceTetQuadrature.hpp"
00036 #include "SundanceBrickQuadrature.hpp"
00037
00038 using namespace Sundance;
00039 using namespace Teuchos;
00040
00041 GaussianQuadrature::GaussianQuadrature(int order)
00042 : QuadratureFamilyBase(order)
00043 {
00044
00045 }
00046
00047 XMLObject GaussianQuadrature::toXML() const
00048 {
00049 XMLObject rtn("GaussianQuadrature");
00050 rtn.addAttribute("order", Teuchos::toString(order()));
00051 return rtn;
00052 }
00053
00054
00055
00056 void GaussianQuadrature::getLineRule(Array<Point>& quadPoints,
00057 Array<double>& quadWeights) const
00058 {
00059 int p = order() + 1;
00060 p = p + (p%2);
00061 int n = p/2;
00062
00063 quadPoints.resize(n);
00064 quadWeights.resize(n);
00065
00066 Gauss1D q1(n, 0.0, 1.0);
00067
00068 for (int i=0; i<n; i++)
00069 {
00070 quadWeights[i] = q1.weights()[i];
00071 quadPoints[i] = Point(q1.nodes()[i]);
00072 }
00073 }
00074
00075 void GaussianQuadrature::getTriangleRule(Array<Point>& quadPoints,
00076 Array<double>& quadWeights) const
00077 {
00078 Array<double> x;
00079 Array<double> y;
00080 Array<double> w;
00081
00082 TriangleQuadrature::getPoints(order(), w, x, y);
00083 quadPoints.resize(w.length());
00084 quadWeights.resize(w.length());
00085 for (int i=0; i<w.length(); i++)
00086 {
00087 quadWeights[i] = 0.5*w[i];
00088 quadPoints[i] = Point(x[i], y[i]);
00089 }
00090 }
00091
00092 void GaussianQuadrature::getQuadRule(Array<Point>& quadPoints,
00093 Array<double>& quadWeights) const
00094 {
00095 Array<double> x;
00096 Array<double> y;
00097 Array<double> w;
00098
00099 QuadQuadrature::getPoints(order(), w, x, y);
00100 quadPoints.resize(w.length());
00101 quadWeights.resize(w.length());
00102 for (int i=0; i<w.length(); i++)
00103 {
00104 quadWeights[i] = w[i];
00105 quadPoints[i] = Point(x[i], y[i]);
00106 }
00107 }
00108
00109 void GaussianQuadrature::getTetRule(Array<Point>& quadPoints,
00110 Array<double>& quadWeights) const
00111 {
00112 Array<double> x;
00113 Array<double> y;
00114 Array<double> z;
00115 Array<double> w;
00116
00117 TetQuadrature::getPoints(order(), w, x, y, z);
00118 quadPoints.resize(w.length());
00119 quadWeights.resize(w.length());
00120 for (int i=0; i<w.length(); i++)
00121 {
00122 quadWeights[i] = w[i]/6.0;
00123 quadPoints[i] = Point(x[i], y[i], z[i]);
00124 }
00125 }
00126
00127 void GaussianQuadrature::getBrickRule(Array<Point>& quadPoints,
00128 Array<double>& quadWeights) const
00129 {
00130 Array<double> x;
00131 Array<double> y;
00132 Array<double> z;
00133 Array<double> w;
00134
00135 BrickQuadrature::getPoints(order(), w, x, y, z);
00136 quadPoints.resize(w.length());
00137 quadWeights.resize(w.length());
00138 for (int i=0; i<w.length(); i++)
00139 {
00140 quadWeights[i] = w[i];
00141 quadPoints[i] = Point(x[i], y[i], z[i]);
00142 }
00143 }
00144
00145 void GaussianQuadrature::getAdaptedWeights(const CellType& cellType ,
00146 int cellDim,
00147 int celLID ,
00148 int facetIndex ,
00149 const Mesh& mesh ,
00150 const ParametrizedCurve& globalCurve ,
00151 Array<Point>& quadPoints ,
00152 Array<double>& quadWeights ,
00153 bool &isCut) const {
00154
00155 isCut = true;
00156 if (mesh.IsSpecialWeightValid() && mesh.hasSpecialWeight(cellDim, celLID))
00157 {
00158 mesh.getSpecialWeight(cellDim, celLID, quadWeights);
00159
00160 return;
00161 }
00162
00163
00164 if (quadPoints.size() <= 1) getPoints(cellType, quadPoints, quadWeights);
00165
00166
00167 isCut = true;
00168 Array<Point> tmpPoints = quadPoints;
00169
00170 Array<int> celLIDs(1);
00171 celLIDs[0] = celLID;
00172
00173 mesh.pushForward( cellDim, celLIDs, quadPoints, tmpPoints );
00174
00175 quadWeights.resize(tmpPoints.size());
00176
00177 for (int i=0 ; i < tmpPoints.size() ; i++){
00178 quadWeights[i] = quadWeights[i] * globalCurve.integrationParameter(tmpPoints[i]);
00179 }
00180
00181
00182 mesh.setSpecialWeight(cellDim, celLID, quadWeights);
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 }