00001
00002
00003
00004
00005
00006
00007
00008 #ifndef SUNDANCEGAUSSLOBATTOQUADRATURE_HPP_
00009 #define SUNDANCEGAUSSLOBATTOQUADRATURE_HPP_
00010
00011 #include "SundanceDefs.hpp"
00012 #include "PlayaTabs.hpp"
00013 #include "SundanceQuadratureFamilyBase.hpp"
00014
00015 namespace Sundance {
00016
00017 class GaussLobattoQuadrature : public QuadratureFamilyBase {
00018
00019 public:
00020
00021
00022 GaussLobattoQuadrature(int order);
00023
00024
00025 virtual ~GaussLobattoQuadrature()
00026 {;}
00027
00028
00029 virtual XMLObject toXML() const;
00030
00031
00032 virtual std::string description() const
00033 {
00034 return "GaussLobattoQuadrature[order=" + Teuchos::toString(order()) + "]";
00035 }
00036
00037
00038 GET_RCP(QuadratureFamilyStub);
00039
00040 protected:
00041
00042 virtual void getLineRule(Array<Point>& quadPointsL,
00043 Array<double>& quadWeights) const;
00044
00045
00046 virtual void getTriangleRule(Array<Point>& quadPoints,
00047 Array<double>& quadWeights) const;
00048
00049
00050 virtual void getQuadRule(Array<Point>& quadPoints,
00051 Array<double>& quadWeights) const;
00052
00053
00054 virtual void getTetRule(Array<Point>& quadPoints,
00055 Array<double>& quadWeights) const ;
00056
00057
00058 virtual void getBrickRule(Array<Point>& quadPoints,
00059 Array<double>& quadWeights) const;
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 virtual void getAdaptedWeights(const CellType& cellType, int cellDim,
00075 int cellLID, int facetIndex, const Mesh& mesh,
00076 const ParametrizedCurve& globalCurve,
00077 Array<Point>& quadPoints, Array<double>& quadWeights,
00078 bool &weightsChanged) const;
00079
00080
00081 virtual void getAdaptedQuadWeights(int cellLID, const Mesh& mesh,
00082 const ParametrizedCurve& globalCurve, Array<Point>& quadPoints,
00083 Array<double>& quadWeights, bool& weightsChanged) const;
00084
00085
00086 virtual void getAdaptedQuadWeights_polygon(int cellLID, const Mesh& mesh,
00087 const ParametrizedCurve& globalCurve, Array<Point>& quadPoints,
00088 Array<double>& quadWeights, bool& weightsChanged) const;
00089
00090
00091 virtual void getAdaptedQuadWeights_surf(int cellLID, const Mesh& mesh,
00092 const ParametrizedCurve& globalCurve, Array<Point>& quadPoints,
00093 Array<double>& quadWeights, bool& weightsChanged) const;
00094
00095 private:
00096
00097
00098 void getTriangleQuadPoints(Array<Point>& pnt , Array<double>& weight ) const;
00099
00100
00101
00102
00103
00104 inline double evalLagrangePol(int i , Array<Point>& pnt , double x) const{
00105 int m = pnt.size();
00106 double p = 1.0;
00107 for (int j = 0 ; (j <= i-1)&&(j<m) ; j++){
00108 p = p *( x - pnt[j][0] )/(pnt[i][0] - pnt[j][0]);
00109 }
00110 for (int j = i+1 ; j < m ; j++){
00111 p = p *( x - pnt[j][0] )/(pnt[i][0] - pnt[j][0]);
00112 }
00113 return p;
00114 }
00115
00116
00117
00118 inline void makeInterpolantQuad( double px, double py, double ofx, double ofy,
00119 int nr1DPoints , int nr2DPoints ,
00120 Array<Point>& linePoints ,
00121 Array<Point>& quadPoints ,
00122 Array<double>& pointWeights ,
00123 Array<double>& weightsOut ,
00124 double areFakt ) const {
00125 int xi,yi;
00126 double xc,yc, intR , valTmp , valBasisX , valBasisY ;
00127
00128
00129 if ( fabs(ofx*ofy) < 1e-14 ) return;
00130
00131
00132
00133
00134
00135 for (int i = 0 ; i < nr2DPoints ; i++){
00136
00137 yi = i % nr1DPoints; xi = i / nr1DPoints; intR = 0.0;
00138
00139 for (int q = 0 ; q < pointWeights.size() ; q++){
00140 xc = px + quadPoints[q][0]*ofx;
00141 yc = py + quadPoints[q][1]*ofy;
00142 valBasisX = evalLagrangePol(xi , linePoints , xc);
00143 valBasisY = evalLagrangePol(yi , linePoints , yc);
00144 valTmp = pointWeights[q] * ( valBasisX * valBasisY);
00145 intR = intR + valTmp;
00146
00147
00148
00149
00150 }
00151 intR = intR * ::fabs(ofx*ofy/areFakt);
00152 weightsOut[i] = weightsOut[i] + intR;
00153 }
00154 }
00155
00156
00157
00158
00159 inline void makeInterpolantBrick(int nr1DPoints , int nr3DPoints ,
00160 Array<Point>& linePoints ,
00161 Array<Point>& quadPoints ,
00162 Array<double>& pointWeights ,
00163 Array<double>& weightsOut ,
00164 double areFakt ) const {
00165 int xi,yi,zi;
00166 double xc,yc,zc, intR , valTmp , valBasisX , valBasisY , valBasisZ;
00167
00168 for (int i = 0 ; i < nr3DPoints ; i++){
00169 zi = i/(nr1DPoints*nr1DPoints); yi = (i%(nr1DPoints*nr1DPoints)) / nr1DPoints;
00170 xi = (i%(nr1DPoints*nr1DPoints)) % nr1DPoints; intR = 0.0;
00171
00172 for (int q = 0 ; q < pointWeights.size() ; q++){
00173 xc = quadPoints[q][0]; yc = quadPoints[q][1]; zc = quadPoints[q][2];
00174 valBasisX = evalLagrangePol(xi , linePoints , xc);
00175 valBasisY = evalLagrangePol(yi , linePoints , yc);
00176 valBasisZ = evalLagrangePol(zi , linePoints , zc);
00177 valTmp = pointWeights[q] * ( valBasisX * valBasisY * valBasisZ );
00178 intR = intR + valTmp;
00179
00180
00181
00182
00183 }
00184 intR = intR * areFakt ;
00185 weightsOut[i] = weightsOut[i] + intR;
00186 }
00187 }
00188
00189
00190
00191
00192 int nrPointin1D_;
00193
00194
00195 int verb_;
00196
00197 static const int quadsEdgesPoints[4][2];
00198
00199
00200 static const int edge3DProjection[12];
00201 };
00202
00203 }
00204
00205 #endif