00001
00002
00003
00004
00005
00006
00007
00008 #include "SundanceTrapesoidQuadrature.hpp"
00009
00010 using namespace Sundance;
00011
00012 TrapesoidQuadrature::TrapesoidQuadrature(int resolution) :
00013 QuadratureFamilyBase(1) , resolution_(resolution) {
00014
00015 resolution_ = (resolution_<2) ? 2 : resolution_;
00016 }
00017
00018 XMLObject TrapesoidQuadrature::toXML() const
00019 {
00020 XMLObject rtn("TrapesoidQuadrature");
00021 rtn.addAttribute("order", Teuchos::toString(order()));
00022 return rtn;
00023 }
00024
00025 void TrapesoidQuadrature::getLineRule(
00026 Array<Point>& quadPoints,
00027 Array<double>& quadWeights) const {
00028
00029 int nrVirInnerPoints = resolution_ - 1;
00030 double innerWeights = 1/((double)nrVirInnerPoints);
00031
00032 int nrPoints = nrVirInnerPoints + 1;
00033 quadPoints.resize(nrPoints);
00034 quadWeights.resize(nrPoints);
00035
00036
00037 quadPoints[0] = Point(0.0);
00038 quadWeights[0] = innerWeights/2.0;
00039
00040 for (int tmp = 1 ; tmp < nrPoints-1 ; tmp++){
00041 quadPoints[tmp] = Point( (double)tmp/(double)(nrPoints-1) );
00042 quadWeights[tmp] = innerWeights;
00043 }
00044
00045
00046 quadPoints[nrPoints-1] = Point(1.0);
00047 quadWeights[nrPoints-1] = innerWeights/2.0;
00048
00049 }
00050
00051
00052 void TrapesoidQuadrature::getTriangleRule(
00053 Array<Point>& quadPoints,
00054 Array<double>& quadWeights) const {
00055
00056 SUNDANCE_ERROR("Triangle rule not available for " << toXML());
00057 }
00058
00059
00060 void TrapesoidQuadrature::getQuadRule(
00061 Array<Point>& quadPoints,
00062 Array<double>& quadWeights) const {
00063
00064 Array<Point> quadPointsLine;
00065 Array<double> quadWeightsLine;
00066
00067 this->getLineRule( quadPointsLine, quadWeightsLine );
00068
00069 int nrPointPerAxis = quadPointsLine.size();
00070
00071 quadPoints.resize(nrPointPerAxis*nrPointPerAxis);
00072 quadWeights.resize(nrPointPerAxis*nrPointPerAxis);
00073
00074 int pcount = 0;
00075 for (int ix = 0 ; ix < nrPointPerAxis ; ix ++){
00076 for (int iy = 0 ; iy < nrPointPerAxis ; iy ++){
00077
00078 quadPoints[pcount] = Point( quadPointsLine[ix][0] , quadPointsLine[iy][0] );
00079 quadWeights[pcount] = quadWeightsLine[ix] * quadWeightsLine[iy];
00080 pcount++;
00081 }
00082 }
00083
00084 }
00085
00086
00087 void TrapesoidQuadrature::getTetRule(
00088 Array<Point>& quadPoints,
00089 Array<double>& quadWeights) const {
00090
00091 SUNDANCE_ERROR("Tet rule not available for " << toXML());
00092 }
00093
00094
00095 void TrapesoidQuadrature::getBrickRule(
00096 Array<Point>& quadPoints,
00097 Array<double>& quadWeights) const {
00098
00099 Array<Point> quadPointsLine;
00100 Array<double> quadWeightsLine;
00101
00102 getLineRule( quadPointsLine, quadWeightsLine );
00103
00104 int nrPointPerAxis = quadPointsLine.size();
00105
00106 quadPoints.resize(nrPointPerAxis*nrPointPerAxis*nrPointPerAxis);
00107 quadWeights.resize(nrPointPerAxis*nrPointPerAxis*nrPointPerAxis);
00108
00109 int pcount = 0;
00110 for (int ix = 0 ; ix < nrPointPerAxis ; ix ++){
00111 for (int iy = 0 ; iy < nrPointPerAxis ; iy ++){
00112 for (int iz = 0 ; iz < nrPointPerAxis ; iz ++){
00113
00114 quadPoints[pcount] = Point( quadPointsLine[ix][0] , quadPointsLine[iy][0] , quadPointsLine[iz][0]);
00115 quadWeights[pcount] = quadWeightsLine[ix] * quadWeightsLine[iy] * quadWeightsLine[iz];
00116 pcount++;
00117 }
00118 }
00119 }
00120 }
00121
00122
00123 void TrapesoidQuadrature::getAdaptedWeights(
00124 const CellType& cellType, int cellDim,
00125 int cellLID, int facetIndex, const Mesh& mesh,
00126 const ParametrizedCurve& globalCurve,
00127 Array<Point>& quadPoints, Array<double>& quadWeights,
00128 bool &weightsChanged) const {
00129
00130 weightsChanged = true;
00131
00132 if (mesh.IsSpecialWeightValid() && mesh.hasSpecialWeight(cellDim, cellLID))
00133 {
00134 mesh.getSpecialWeight(cellDim, cellLID, quadWeights);
00135
00136 return;
00137 }
00138
00139
00140 if (quadPoints.size() <= 1) getPoints(cellType, quadPoints, quadWeights);
00141
00142
00143 weightsChanged = true;
00144 Array<Point> tmpPoints = quadPoints;
00145
00146 Array<int> celLIDs(1);
00147 celLIDs[0] = cellLID;
00148
00149 mesh.pushForward( cellDim, celLIDs, quadPoints, tmpPoints );
00150
00151 quadWeights.resize(tmpPoints.size());
00152
00153 for (int i=0 ; i < tmpPoints.size() ; i++){
00154 quadWeights[i] = quadWeights[i] * globalCurve.integrationParameter(tmpPoints[i]);
00155 }
00156
00157
00158 mesh.setSpecialWeight(cellDim, cellLID, quadWeights);
00159 }