Sundance 2.1

00001 /*
00002  * SundanceGaussLobattoQuadrature.cpp
00003  *
00004  *  Created on: Jan 20, 2011
00005  *      Author: benk
00006  */
00007 
00008 #include "SundanceGaussLobattoQuadrature.hpp"
00009 #include "SundanceCurveIntegralCalc.hpp"
00010 #include "SundancePolygon2D.hpp"
00011 #include "SundanceSurf3DCalc.hpp"
00012 #include "PlayaTabs.hpp"
00013 
00014 using namespace Sundance;
00015 
00016 
00017 int const GaussLobattoQuadrature::quadsEdgesPoints[4][2] = { {0,1} , {0,2} , {1,3} , {2,3} };
00018                                                        //  0   1   2   3   4   5   6   7   8   9  10  11
00019 int const GaussLobattoQuadrature::edge3DProjection[12] = { 0 , 1 , 2 , 1 , 2 , 0 , 2 , 2 , 0 , 1 , 1 , 0};
00020 
00021 
00022 GaussLobattoQuadrature::GaussLobattoQuadrature(int order) :
00023    QuadratureFamilyBase(order) {
00024   nrPointin1D_ = order+1;
00025   verb_ = 0;
00026 }
00027 
00028 XMLObject GaussLobattoQuadrature::toXML() const
00029 {
00030   XMLObject rtn("GaussLobattoQuadrature");
00031   rtn.addAttribute("order", Teuchos::toString(order()));
00032   return rtn;
00033 }
00034 
00035 
00036 void GaussLobattoQuadrature::getTriangleRule(
00037     Array<Point>& quadPoints,
00038     Array<double>& quadWeights) const {
00039     // todo: implement this
00040   SUNDANCE_ERROR("Triangle rule not available for " << toXML());
00041 }
00042 
00043 void GaussLobattoQuadrature::getQuadRule(
00044     Array<Point>& quadPoints,
00045     Array<double>& quadWeights) const {
00046   Array<Point> quadPointsLine;
00047   Array<double> quadWeightsLine;
00048   // get the line rule
00049     this->getLineRule( quadPointsLine, quadWeightsLine );
00050 
00051     int nrPointPerAxis = quadPointsLine.size();
00052     // we simply take the tensor product
00053     quadPoints.resize(nrPointPerAxis*nrPointPerAxis);
00054     quadWeights.resize(nrPointPerAxis*nrPointPerAxis);
00055 
00056     int pcount = 0;
00057     for (int ix = 0 ; ix < nrPointPerAxis ; ix ++){
00058       for (int iy = 0 ; iy < nrPointPerAxis ; iy ++){
00059         // here we take the tensor product of the
00060         quadPoints[pcount] = Point( quadPointsLine[ix][0] , quadPointsLine[iy][0] );
00061         quadWeights[pcount] = quadWeightsLine[ix] * quadWeightsLine[iy];
00062         pcount++;
00063       }
00064     }
00065 }
00066 
00067 
00068 void GaussLobattoQuadrature::getTetRule(
00069     Array<Point>& quadPoints,
00070     Array<double>& quadWeights) const {
00071    // todo: implement this
00072    SUNDANCE_ERROR("Tet rule not available for " << toXML());
00073 }
00074 
00075 
00076 void GaussLobattoQuadrature::getBrickRule(
00077     Array<Point>& quadPoints,
00078     Array<double>& quadWeights) const {
00079   Array<Point> quadPointsLine;
00080   Array<double> quadWeightsLine;
00081   // get the line rule
00082     this->getLineRule( quadPointsLine, quadWeightsLine );
00083 
00084     int nrPointPerAxis = quadPointsLine.size();
00085     // we simply take the tensor product
00086     quadPoints.resize(nrPointPerAxis*nrPointPerAxis*nrPointPerAxis);
00087     quadWeights.resize(nrPointPerAxis*nrPointPerAxis*nrPointPerAxis);
00088 
00089     int pcount = 0;
00090     for (int iz = 0 ; iz < nrPointPerAxis ; iz ++){
00091     for (int iy = 0 ; iy < nrPointPerAxis ; iy ++){
00092         for (int ix = 0 ; ix < nrPointPerAxis ; ix ++){
00093         // here we take the tensor product of the
00094         quadPoints[pcount] = Point( quadPointsLine[ix][0] , quadPointsLine[iy][0] , quadPointsLine[iz][0]);
00095         quadWeights[pcount] = quadWeightsLine[ix] * quadWeightsLine[iy] * quadWeightsLine[iz];
00096         pcount++;
00097       }
00098       }
00099     }
00100 }
00101 
00102 
00103 void GaussLobattoQuadrature::getAdaptedWeights(
00104     const CellType& cellType, int cellDim,
00105     int cellLID, int facetIndex, const Mesh& mesh,
00106     const ParametrizedCurve& globalCurve,
00107     Array<Point>& quadPoints, Array<double>& quadWeights,
00108     bool &weightsChanged) const {
00109 
00110   Tabs tabs;
00111 
00112   weightsChanged = true; // todo: this not always might be true to save computations we might test if this is the case
00113   if (mesh.IsSpecialWeightValid() && mesh.hasSpecialWeight(cellDim, cellLID))
00114   {
00115     mesh.getSpecialWeight(cellDim, cellLID, quadWeights);
00116     //SUNDANCE_MSG3(verb, tabs << "GaussianQuadrature::getAdaptedWeights Found cached weights for cell LID " << cellLID)
00117     return;
00118   }
00119   // if we have no quad points then get them
00120   if (quadPoints.size() <= 1) getPoints(cellType,  quadPoints, quadWeights);
00121 
00122   // Maximal cell type
00123   CellType maxCellType = mesh.cellType(mesh.spatialDim());
00124 
00125   switch (maxCellType)
00126   {
00127   // We have a triangle based mesh
00128   case TriangleCell:
00129   {
00130     switch (cellType)
00131     {
00132     case TriangleCell:
00133     {
00134       // todo: implement this
00135       break;
00136     }
00137     default:
00138 #ifndef TRILINOS_7
00139       SUNDANCE_ERROR ("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh")
00140       ;
00141 #else
00142       SUNDANCE_ERROR7("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh");
00143 #endif
00144     }
00145     break;
00146   }
00147   case QuadCell:
00148   {
00149     switch(cellType)
00150     {
00151     case QuadCell:
00152     {
00153       // call the method to integrate the Quad
00154 
00155       if (CurveIntegralCalc::usePolygoneCurveQuadrature( cellType , mesh , globalCurve)){
00156         // if the curve is a polygon then use a different
00157         getAdaptedQuadWeights_polygon(cellLID, mesh, globalCurve, quadPoints, quadWeights, weightsChanged);
00158       }
00159       else
00160       {
00161         getAdaptedQuadWeights(cellLID, mesh, globalCurve, quadPoints, quadWeights, weightsChanged);
00162       }
00163       break;
00164     }
00165     default:
00166 #ifndef TRILINOS_7
00167     SUNDANCE_ERROR ("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh")
00168     ;
00169 #else
00170     SUNDANCE_ERROR7("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh");
00171 #endif
00172     }
00173     break;
00174   }
00175   case BrickCell:
00176   {
00177     switch(cellType)
00178     {
00179     case BrickCell:
00180     {
00181       // call the method which does the job for 3D
00182       getAdaptedQuadWeights_surf(cellLID, mesh, globalCurve, quadPoints, quadWeights, weightsChanged);
00183       break;
00184     }
00185     default: {}
00186       }
00187     break;
00188   }
00189   default:
00190 #ifndef TRILINOS_7
00191     SUNDANCE_ERROR("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for cell type " << cellType)
00192     ;
00193 #else
00194     SUNDANCE_ERROR7("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for cell type " << cellType);
00195 #endif
00196   }
00197 
00198   //store the weights in the mesh
00199   mesh.setSpecialWeight(cellDim, cellLID, quadWeights);
00200 }
00201 
00202 
00203 
00204 void GaussLobattoQuadrature::getAdaptedQuadWeights_polygon(int cellLID, const Mesh& mesh,
00205     const ParametrizedCurve& globalCurve, Array<Point>& quadPoints,
00206     Array<double>& quadWeights, bool& weightsChanged) const {
00207 
00208   // this works only in 2D for polygons, it is important that we ask the underlying object and not the object itself
00209   const Polygon2D* poly = dynamic_cast<const Polygon2D* > ( globalCurve.ptr().get()->getUnderlyingCurveObj() );
00210   if ( poly == 0 ) {
00211     TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error , " getAdaptedQuadWeights_polygon , paramCurve must be a Polygon and the cell must be a quad" );
00212     return;
00213   }
00214 
00215   // the main idea of the algorithm is to get a raw of points
00216   // (P1,P2,P3,P4,P5) and (alpha1,aplpha2), the coefficients
00217   // this raw of points must be ordered correctly and contains intersection and polygon points
00218   // by integrating the left side of this curve we can calculate the exact integral
00219 
00220   // - get the intersection points
00221   // - get the polygon points
00222   // form out of these points a row of points with the shortest distance
00223   // choose as starting point the one with the smallest x and y coordinate
00224 
00225   // use this raw of points for quadrature
00226   //  - see in which direction is always increasing the coordinate
00227   //  - choose one direction (x or y) to parse the raw of points (if not always increasing then take the sign in consideration)
00228 
00229   Array<Point> quadOfPoints(4);
00230   Array<Point> intesectPoints(8);
00231   Array<Point> pointStack(0);
00232   Array<bool>  isInnerPoint(0);
00233   Array<Point> pointStackTmp(0);
00234   Array<Point> allPolygonPoints(0);
00235   int nrIPoints , pointsNr , polyNrPoint ;
00236   double eps = 1e-14;
00237 
00238   SUNDANCE_MSG3(verb_, " ---------- START METHOD -----------  ");
00239 
00240   // get the points of the quad
00241   for (int i = 0 ; i < 4 ; i++){
00242     quadOfPoints[i] = mesh.nodePosition( mesh.facetLID(2,cellLID , 0 , i , nrIPoints ) );
00243   }
00244 
00245   // get all the intersection points
00246   pointsNr = 0;
00247   for (int i = 0 ; i < 4 ; i++)
00248   {
00249     // for each edge test intersection point
00250     globalCurve.returnIntersectPoints( quadOfPoints[quadsEdgesPoints[i][0]] , quadOfPoints[quadsEdgesPoints[i][1]],
00251          nrIPoints , intesectPoints );
00252     //SUNDANCE_MSG3(verb_, " Test edge: "<< quadOfPoints[quadsEdgesPoints[i][0]] << " -- " << quadOfPoints[quadsEdgesPoints[i][1]] );
00253     for (int j = 0 ; j < nrIPoints ; j++){
00254       bool addThisPoint = true;
00255       // add these points only when they are not duplicate
00256       for (int tmpP = 0 ; tmpP < pointsNr ; tmpP++){
00257         double dist = ::sqrt( (pointStackTmp[tmpP] - intesectPoints[j])*(pointStackTmp[tmpP] - intesectPoints[j]) );
00258         if ( dist < eps) addThisPoint = false;
00259       }
00260       // if the point is not a duplicate then add
00261       if (addThisPoint) {
00262         pointStack.resize( pointsNr + 1  );
00263         pointStackTmp.resize( pointsNr + 1 );
00264         isInnerPoint.resize( pointsNr + 1 );
00265         // transform to reference coordinates back
00266         pointStack[pointsNr] = Point( ( intesectPoints[j][0] - quadOfPoints[0][0])/( quadOfPoints[3][0]-quadOfPoints[0][0] ) ,
00267                                 ( intesectPoints[j][1] - quadOfPoints[0][1])/( quadOfPoints[3][1]-quadOfPoints[0][1] )  );
00268         // we store the intersection point also in real coordinates
00269         pointStackTmp[pointsNr] = intesectPoints[j];
00270         isInnerPoint[pointsNr] = false;
00271         //SUNDANCE_MSG3(verb_, " adding Ints Points:"<< pointsNr << " P:"
00272         //    << pointStack[pointsNr] << " Real P:" << intesectPoints[j] );
00273         pointsNr = pointsNr + 1;
00274       }
00275     }
00276   }
00277 
00278   // get the point
00279   poly->getCellsPolygonesPoint( cellLID ,allPolygonPoints);
00280   polyNrPoint = allPolygonPoints.size();
00281   //SUNDANCE_MSG3(verb_, " Polygon has :"<< polyNrPoint << " points inside the cell");
00282   // add also the points which come from the polygon
00283   for (int j = 0 ; j < polyNrPoint ; j++){
00284     bool addThisPoint = true;
00285     // add these points only when they are not duplicate
00286     for (int tmpP = 0 ; tmpP < pointsNr ; tmpP++){
00287       double dist = ::sqrt( (pointStackTmp[tmpP] - allPolygonPoints[j])*(pointStackTmp[tmpP] - allPolygonPoints[j]) );
00288       if ( dist < eps) addThisPoint = false;
00289     }
00290     // if the point is not a duplicate then add
00291     if (addThisPoint) {
00292       pointStack.resize( pointsNr + 1  );
00293       pointStackTmp.resize( pointsNr + 1 );
00294       isInnerPoint.resize( pointsNr + 1 );
00295       // transform to reference coordinates back
00296       pointStack[pointsNr] = Point( ( allPolygonPoints[j][0] - quadOfPoints[0][0])/( quadOfPoints[3][0]-quadOfPoints[0][0] ) ,
00297                               ( allPolygonPoints[j][1] - quadOfPoints[0][1])/( quadOfPoints[3][1]-quadOfPoints[0][1] )  );
00298       // we store the intersection point also in real coordinates
00299       pointStackTmp[pointsNr] = allPolygonPoints[j];
00300       isInnerPoint[pointsNr] = true;
00301       //SUNDANCE_MSG3(verb_, " adding Polygon Points:"<< pointsNr << " P:"
00302       //    << pointStack[pointsNr] << " Real P:" << allPolygonPoints[j] );
00303       pointsNr = pointsNr + 1;
00304     }
00305   }
00306 
00307   // now we have all the points we just need to get them in a sorted way, after distance (Manhattan distance)
00308   // first we select the point which is closes to the (0,0) Point and is an intersection point of the edges
00309   double dist = 1e+100 , firstind = -1;
00310   for (int tmpP = 0 ; tmpP < pointsNr ; tmpP++){
00311     // we look only for intersection points to have as starting point
00312     if ( (pointStack[tmpP][0] + pointStack[tmpP][1] < dist) && (!isInnerPoint[tmpP])){
00313       dist = pointStack[tmpP][0] + pointStack[tmpP][1];
00314       firstind = tmpP;
00315     }
00316   }
00317 
00318   Array<Point> orderedPoints(1);
00319   Array<Point> orderedPointsRealCoords(1);
00320   Array<bool> orderedPointsIsIntersection(1);
00321   orderedPoints[0] = pointStack[firstind];
00322   orderedPointsRealCoords[0] = pointStackTmp[firstind];
00323   orderedPointsIsIntersection[0] = (!isInnerPoint[firstind]);
00324   //SUNDANCE_MSG3(verb_, "Add first point to position 0  from firstind:"<< firstind << " , P:" << orderedPoints[0]);
00325 
00326   int orderedP = 1;
00327 
00328   // with this for loop
00329   for (int tmpP = 0 ; tmpP < pointsNr ; tmpP++){
00330     // the first point which we added we do not want to add again
00331     if (tmpP == firstind) continue;
00332 
00333     // add all the 1..pointsNr points to the ordered list
00334     double dist = 1e+100 , dist_tmp;
00335     int ind = -1;
00336     //SUNDANCE_MSG3(verb_, " Trying to insert: "<< tmpP << " out of " << pointsNr);
00337     for (int p = 0 ; p <= orderedP ; p++){
00338       dist_tmp = dist;
00339       //SUNDANCE_MSG3(verb_, " probe position " << p << " , orderedP:"<<orderedP);
00340       if (( p == orderedP ) || (p == 0) ){
00341         if (!isInnerPoint[tmpP]){
00342           Point poin = orderedPoints[(p == orderedP)?p-1:p] - pointStack[tmpP];
00343           // check this, what kind of distance * 1.0 should be here added
00344           dist_tmp = 1.0 * ::sqrt(poin*poin);
00345         }
00346       }
00347       else {
00348         Point poin1 = orderedPoints[p-1] - pointStack[tmpP];
00349         Point poin2 = orderedPoints[p] - pointStack[tmpP];
00350         Point prev_dist = orderedPoints[p] - orderedPoints[p-1];
00351         dist_tmp = ::sqrt(poin1*poin1) + ::sqrt(poin2*poin2) - ::sqrt(prev_dist*prev_dist);
00352       }
00353       //SUNDANCE_MSG3(verb_, " dist_tmp: " << dist_tmp << " , dist:" << dist );
00354       // if this distance is shorter than mark this position
00355       if (dist_tmp < dist){
00356         dist = dist_tmp;
00357         ind = p;
00358       }
00359     }
00360 
00361     // the first point should be inserted always to the back
00362     if ( (ind == 0) && ( orderedP <= 1) ){ ind = 1;}
00363 
00364     // if "ind"
00365     TEUCHOS_TEST_FOR_EXCEPTION( ind < 0 , std::runtime_error , " GaussLobattoQuadrature::getAdaptedQuadWeights_polygon ind < 0 ind:"<<ind );
00366     // insert the point to the "ind" position
00367     orderedPoints.resize(orderedP+1);
00368     orderedPointsRealCoords.resize(orderedP+1);
00369     orderedPointsIsIntersection.resize(orderedP+1);
00370     for (int p = orderedP-1 ; p >= ind ; p--){
00371       orderedPoints[p+1] = orderedPoints[p];
00372       orderedPointsRealCoords[p+1] = orderedPointsRealCoords[p];
00373       orderedPointsIsIntersection[p+1] = orderedPointsIsIntersection[p];
00374     }
00375     //SUNDANCE_MSG3(verb_, "Add point to position ind:"<< ind <<  " , tmpP:" << tmpP << " , orderedP:" << orderedP );
00376     //SUNDANCE_MSG3(verb_, "pointStack[tmpP]:"<< pointStack[tmpP] <<  " , pointStackTmp[tmpP]:" << pointStackTmp[tmpP] );
00377     orderedPoints[ind] = pointStack[tmpP];
00378     orderedPointsRealCoords[ind] = pointStackTmp[tmpP];
00379     orderedPointsIsIntersection[ind] = (!isInnerPoint[tmpP]);
00380     orderedP = orderedP + 1;
00381   }
00382 
00383   // print the points which were ordered
00384   for (int pi = 0 ; pi < orderedP ; pi++){
00385     //SUNDANCE_MSG3(verb_, "Ordered Points pi:"<< pi << " P:" << orderedPoints[pi] );
00386   }
00387   // at this stage we have the ordered raw of point we just have to use it
00388   // the ordered
00389   bool increaseX = true;
00390   bool increaseY = true;
00391   double minX = orderedPoints[0][0], maxX = orderedPoints[0][0], minY = orderedPoints[0][1], maxY = orderedPoints[0][1];
00392   for (int p = 1 ; p < orderedP ; p++){
00393     if ( orderedPoints[p-1][0] - eps > orderedPoints[p][0]) increaseX = false;
00394     if ( orderedPoints[p-1][1] - eps > orderedPoints[p][1]) increaseY = false;
00395     minX = (orderedPoints[p][0] < minX) ? orderedPoints[p][0] : minX;
00396     maxX = (orderedPoints[p][0] > maxX) ? orderedPoints[p][0] : maxX;
00397     minY = (orderedPoints[p][1] < minY) ? orderedPoints[p][1] : minY;
00398     maxY = (orderedPoints[p][1] > maxY) ? orderedPoints[p][1] : maxY;
00399   }
00400 
00401   // get all the possible weights and quadrature points
00402   Array<Point> linePoints;
00403   Array<double> lineWeights;
00404   getLineRule( linePoints , lineWeights);
00405   int nr1DPoints = linePoints.size();
00406   int nr2DPoints = nr1DPoints * nr1DPoints;
00407   Array<Point> quadQuadPoints;
00408   Array<double> quadQuadWeights;
00409   getQuadRule( quadQuadPoints , quadQuadWeights);
00410   Array<Point> trianglePoints;
00411   Array<double> triangleWeights;
00412   getTriangleQuadPoints( trianglePoints , triangleWeights);
00413 
00414   // the integration coefficients , we get the coefficient left from the curve and right from the curve
00415   double alphaLeft , alphaRight;
00416   Point curvderi = orderedPointsRealCoords[0] - orderedPointsRealCoords[1];
00417   Point norm_vec = Point( -curvderi[1] , curvderi[0] );
00418   norm_vec = norm_vec / ::sqrt(norm_vec*norm_vec);
00419   double relative_length = 1e-2 * ::sqrt(curvderi*curvderi);
00420   //SUNDANCE_MSG3(verb_, " First segment P0: " << orderedPointsRealCoords[0] << " , P1:" <<orderedPointsRealCoords[1]);
00421   //SUNDANCE_MSG3(verb_, " Alpha left P: " << (0.5*orderedPointsRealCoords[0] + 0.5*orderedPointsRealCoords[1]) + relative_length*norm_vec );
00422   //SUNDANCE_MSG3(verb_, " Alpha right P: " << (0.5*orderedPointsRealCoords[0] + 0.5*orderedPointsRealCoords[1]) - relative_length*norm_vec );
00423   alphaLeft = globalCurve.integrationParameter( (0.5*orderedPointsRealCoords[0] + 0.5*orderedPointsRealCoords[1]) + relative_length*norm_vec );
00424   alphaRight = globalCurve.integrationParameter( (0.5*orderedPointsRealCoords[0] + 0.5*orderedPointsRealCoords[1]) - relative_length*norm_vec );
00425 
00426   //SUNDANCE_MSG3(verb_, "minX:"<< minX << " , maxX:" << maxX << " , minY:" << minY << " , maxY:" << maxY <<
00427   //    " , alphaLeft:" << alphaLeft << " , alphaRight:" << alphaRight);
00428   //SUNDANCE_MSG3(verb_, "increaseX:"<< increaseX << " , increaseY:" << increaseY );
00429 
00430   // first quadrate the whole quad
00431   Array<double> wholeWeightsQuad( quadWeights.size() , 0.0 );
00432     makeInterpolantQuad( 0.0, 0.0 , 1.0 , 1.0, nr1DPoints , nr2DPoints ,linePoints ,
00433                      quadQuadPoints , quadQuadWeights ,wholeWeightsQuad , 1.0 );
00434 
00435   Array<double> leftWeightsQuad( quadWeights.size() , 0.0 );
00436   // here we have a big branch, eighter to parse in the X direction or in the Y direction
00437   if ( (increaseX) || (increaseY == false)){
00438     // quadrate in the X direction
00439     // first make a quadrature from 0 to MinX
00440     //SUNDANCE_MSG3(verb_, " Increase X , integrate initial quad");
00441       makeInterpolantQuad( 0.0, orderedPoints[0][1] ,
00442                        minX , 1.0-orderedPoints[0][1],
00443                        nr1DPoints , nr2DPoints ,linePoints ,
00444                        quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00445 
00446       // for each segment make a triangle and quad quadrature
00447       for (int seg = 1 ; seg < orderedP ; seg++){
00448         //SUNDANCE_MSG3(verb_, " Integrate segment : " << seg );
00449         Point p0 = orderedPoints[seg-1];
00450         Point p1 = orderedPoints[seg];
00451         if (p0[0] < p1[0] ){
00452           // X is increasing
00453           if (p0[1] < p1[1]){
00454             // Y increasing
00455               makeInterpolantQuad( p0[0] , p1[1] ,
00456                                p1[0] - p0[0] , 1.0-p1[1] ,
00457                                nr1DPoints , nr2DPoints ,linePoints ,
00458                                quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00459               makeInterpolantQuad( p0[0] , p1[1] ,
00460                                    p1[0] - p0[0] , p0[1] - p1[1] ,
00461                                nr1DPoints , nr2DPoints ,linePoints ,
00462                                trianglePoints , triangleWeights ,leftWeightsQuad , 2.0  );
00463           }else{
00464             // Y decreasing
00465               makeInterpolantQuad( p0[0] , p0[1] ,
00466                                p1[0] - p0[0] , 1.0-p0[1] ,
00467                                nr1DPoints , nr2DPoints ,linePoints ,
00468                                quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00469               makeInterpolantQuad( p1[0] , p0[1] ,
00470                                    p0[0] - p1[0] , p1[1] - p0[1] ,
00471                                nr1DPoints , nr2DPoints ,linePoints ,
00472                                trianglePoints , triangleWeights ,leftWeightsQuad , 2.0  );
00473           }
00474         }
00475         else{
00476           // X is decreasing (this should happen rather seldom)
00477           if (p0[1] < p1[1]){
00478             // Y increasing
00479               makeInterpolantQuad( p1[0] , 0.0 ,
00480                                p1[0] - p0[0] , p0[1] ,
00481                                nr1DPoints , nr2DPoints ,linePoints ,
00482                                quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00483               makeInterpolantQuad( p1[0] , p0[1] ,
00484                                    p0[0] - p1[0] , p1[1] - p0[1] ,
00485                                nr1DPoints , nr2DPoints ,linePoints ,
00486                                trianglePoints , triangleWeights ,leftWeightsQuad , 2.0  );
00487           }else{
00488             // Y decreasing
00489               makeInterpolantQuad( p1[0] , 0.0 ,
00490                                p0[0] - p1[0] , p1[1] ,
00491                                nr1DPoints , nr2DPoints ,linePoints ,
00492                                quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00493               makeInterpolantQuad( p0[0] , p1[1] ,
00494                                    p1[0] - p0[0] , p0[1] - p1[1] ,
00495                                nr1DPoints , nr2DPoints ,linePoints ,
00496                                trianglePoints , triangleWeights ,leftWeightsQuad , 2.0  );
00497           }
00498         }
00499       }
00500 
00501       //SUNDANCE_MSG3(verb_, " Integrate rest of the quad ");
00502       // quadrate the rest of the quad
00503       makeInterpolantQuad( maxX , orderedPoints[orderedP-1][1] ,
00504                        1.0-maxX , 1.0-orderedPoints[orderedP-1][1] ,
00505                            nr1DPoints , nr2DPoints ,linePoints ,
00506                        quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00507   }
00508   else{
00509     //SUNDANCE_MSG3(verb_, " Increase Y , integrate initial quad");
00510     // quadrate in Y direction
00511     // first make a quadrature from 0 to MinY
00512       makeInterpolantQuad( 0.0 , 0.0 ,
00513                        orderedPoints[0][0] , minY,
00514                        nr1DPoints , nr2DPoints ,linePoints ,
00515                        quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00516 
00517       // for each segment make a triangle and quad quadrature
00518       for (int seg = 1 ; seg < orderedP ; seg++){
00519         //SUNDANCE_MSG3(verb_, " Integrate segment : " << seg );
00520         Point p0 = orderedPoints[seg-1];
00521         Point p1 = orderedPoints[seg];
00522         // Y must be always increasing !!!
00523         if (p0[0] < p1[0]){
00524           // X increasing
00525             makeInterpolantQuad( 0.0 , p0[1] ,
00526                              p0[0] , p1[1]-p0[1] ,
00527                              nr1DPoints , nr2DPoints ,linePoints ,
00528                              quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00529             makeInterpolantQuad( p0[0] , p1[1] ,
00530                                  p1[0] - p0[0] , p0[1] - p1[1] ,
00531                              nr1DPoints , nr2DPoints ,linePoints ,
00532                              trianglePoints , triangleWeights ,leftWeightsQuad , 2.0  );
00533         }else{
00534           // X decreasing
00535             makeInterpolantQuad( 0.0 , p0[1] ,
00536                              p1[0] , p1[1]-p0[1] ,
00537                              nr1DPoints , nr2DPoints ,linePoints ,
00538                              quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00539             makeInterpolantQuad( p1[0] , p0[1] ,
00540                                  p0[0] - p1[0] , p1[1] - p0[1] ,
00541                              nr1DPoints , nr2DPoints ,linePoints ,
00542                              trianglePoints , triangleWeights ,leftWeightsQuad , 2.0  );
00543         }
00544       }
00545 
00546       //SUNDANCE_MSG3(verb_, " Integrate rest of the quad ");
00547       // quadrate the rest of the quad
00548       makeInterpolantQuad( 0.0 , maxY ,
00549                        orderedPoints[orderedP-1][0] , 1.0-maxY ,
00550                            nr1DPoints , nr2DPoints ,linePoints ,
00551                        quadQuadPoints , quadQuadWeights ,leftWeightsQuad , 1.0 );
00552   }
00553 
00554 
00555   // here just add up   W(i) = alphaLeft*leftW(i) + alphaRight*(wholeW(i)-leftW(i))
00556   double sumWeights = 0.0;
00557   for (int q = 0 ; q < quadPoints.size() ; q++ ){
00558     quadWeights[q] = alphaLeft * leftWeightsQuad[q] + alphaRight*( wholeWeightsQuad[q]- leftWeightsQuad[q] );
00559     sumWeights = sumWeights + quadWeights[q];
00560   }
00561   SUNDANCE_MSG3(verb_, " ---------- END METHOD -----------  sum weights = " << sumWeights);
00562 }
00563 
00564 
00565 
00566 void GaussLobattoQuadrature::getAdaptedQuadWeights(int cellLID, const Mesh& mesh,
00567     const ParametrizedCurve& globalCurve, Array<Point>& quadPoints,
00568     Array<double>& quadWeights, bool& weightsChanged) const{
00569 
00570   int maxStack = 1000;
00571   //int maxLevel = 5;
00572   int stackIndex = 0;
00573   Tabs tabs;
00574 
00575   Array< Array<Point> > pointStack(maxStack);     // -> the intersection points of each quad cell in the stack (in ref coords)
00576   Array< Array<Point> > quadPointStack(maxStack); // -> the points which form the quad cell (in ref coords)
00577   Array< Array<int> > intersectionPointStack(maxStack); // -> which edges are intersected for the actual quad cell
00578   Array<int> levelStack(maxStack);                      // -> the refinement level of the actual quad cell
00579   Array<int> refinedStack(maxStack,-1);                 // -> if a quad cell has been replaced then this will be > 0 , -1 otherwise
00580   Array<int> intersectioncaseStack(maxStack,-1);        // -> the intersection case can be from 0,..,6 , if it is something else then refine
00581   Array<double> alpha1(maxStack);                       // -> the integration coefficient for the first part of the integral
00582   Array<double> alpha2(maxStack);                       // -> the integration coefficient for the second part of the integral
00583   // the actual points of the quad cell and the actual intersection points in real coords
00584   Array<Point> quadOfPoints(4);
00585   Array<Point> intesectPoints(8);
00586   int nrIPoints , tmp;
00587 
00588   // first add the initial quad to the quad stack
00589   for (int i = 0 ; i < 4 ; i++){
00590     quadOfPoints[i] = mesh.nodePosition( mesh.facetLID(2,cellLID , 0 , i , nrIPoints ) );
00591   }
00592 
00593   // first divide the parent cells in smaller quads ,
00594   pointStack[0].resize(0);
00595   quadPointStack[0].resize(4);
00596   quadPointStack[0][0] = Point(0.0,0.0); quadPointStack[0][1] = Point(1.0,0.0);
00597   quadPointStack[0][2] = Point(0.0,1.0); quadPointStack[0][3] = Point(1.0,1.0);
00598   intersectionPointStack[0].resize(0);
00599   levelStack[0] = 0;
00600   tmp = 0;
00601   // for each edge test for intersection
00602   for (int i = 0 ; i < 4 ; i++)
00603   {
00604     // for each edge test intersection point
00605     globalCurve.returnIntersectPoints( quadOfPoints[quadsEdgesPoints[i][0]] , quadOfPoints[quadsEdgesPoints[i][1]],
00606          nrIPoints , intesectPoints );
00607     SUNDANCE_MSG3(verb_, " Test edge: "<< quadOfPoints[quadsEdgesPoints[i][0]] << " -- " << quadOfPoints[quadsEdgesPoints[i][1]] );
00608     pointStack[0].resize( tmp + nrIPoints );
00609     intersectionPointStack[0].resize( tmp + nrIPoints );
00610     for (int j = 0 ; j < nrIPoints ; j++){
00611       // transform to reference coordinates back
00612       pointStack[0][tmp+j] = Point( ( intesectPoints[j][0] - quadOfPoints[0][0])/( quadOfPoints[3][0]-quadOfPoints[0][0] ) ,
00613                                 ( intesectPoints[j][1] - quadOfPoints[0][1])/( quadOfPoints[3][1]-quadOfPoints[0][1] )  );
00614       intersectionPointStack[0][tmp+j] = i;
00615       SUNDANCE_MSG3(verb_, " adding Ints Points:"<< tmp+j << " edge:" << intersectionPointStack[0][tmp+j] << " P:"
00616           << pointStack[0][tmp+j] << " Real P:" << intesectPoints[j] );
00617     }
00618     tmp = tmp + nrIPoints;
00619   }
00620 
00621   // select the correct case
00622   intersectioncaseStack[0] = -1;
00623   if (intersectionPointStack[0].size() == 0){
00624     intersectioncaseStack[0] = 6; // no intersection point
00625     alpha1[0] = alpha2[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00626   }
00627   if (intersectionPointStack[0].size() == 2){
00628        if ((intersectionPointStack[0][0] == 0 ) && (intersectionPointStack[0][1] == 1 )) {
00629          intersectioncaseStack[0] = 0;
00630          alpha1[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00631          alpha2[0] = globalCurve.integrationParameter(quadOfPoints[3]);
00632        }
00633        if ((intersectionPointStack[0][0] == 0 ) && (intersectionPointStack[0][1] == 2 )) {
00634          intersectioncaseStack[0] = 2;
00635          alpha1[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00636          alpha2[0] = globalCurve.integrationParameter(quadOfPoints[1]);
00637        }
00638        if ((intersectionPointStack[0][0] == 0 ) && (intersectionPointStack[0][1] == 3 )) {
00639          alpha1[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00640          alpha2[0] = globalCurve.integrationParameter(quadOfPoints[3]);
00641          if (pointStack[0][0][0] > pointStack[0][1][0]) intersectioncaseStack[0] = 41;
00642          else intersectioncaseStack[0] = 42;
00643        }
00644        if ((intersectionPointStack[0][0] == 1 ) && (intersectionPointStack[0][1] == 2 )) {
00645          alpha1[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00646          alpha2[0] = globalCurve.integrationParameter(quadOfPoints[3]);
00647          if (pointStack[0][0][1] > pointStack[0][1][1]) intersectioncaseStack[0] = 51;
00648          else intersectioncaseStack[0] = 52;
00649        }
00650        if ((intersectionPointStack[0][0] == 1 ) && (intersectionPointStack[0][1] == 3 )) {
00651          alpha1[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00652          alpha2[0] = globalCurve.integrationParameter(quadOfPoints[2]);
00653          intersectioncaseStack[0] = 1;
00654        }
00655        if ((intersectionPointStack[0][0] == 2 ) && (intersectionPointStack[0][1] == 3 )) {
00656          alpha1[0] = globalCurve.integrationParameter(quadOfPoints[0]);
00657          alpha2[0] = globalCurve.integrationParameter(quadOfPoints[3]);
00658          intersectioncaseStack[0] = 3;
00659        }
00660   }
00661   SUNDANCE_MSG3(verb_, " intersectioncaseStack[0]:"<< intersectioncaseStack[0] << " , alpha1[0]:" << alpha1[0]
00662                          << " , alpha2[0]" << alpha2[0]);
00663   stackIndex = 1;
00664   //quadsEdgesPoints
00665 
00666 
00667   // the criterion for division is if the are after the refinement is different significantly
00668   // then before the division
00669   // put these quads with the intersection points in one list
00670   // todo: here we make no while loop since then also the curve integrals should be adaptive
00671 
00672 
00673   // get all the possible weights and quadrature points
00674   Array<Point> linePoints;
00675   Array<double> lineWeights;
00676   getLineRule( linePoints , lineWeights);
00677   Array<Point> quadQuadPoints;
00678   Array<double> quadQuadWeights;
00679   getQuadRule( quadQuadPoints , quadQuadWeights);
00680   Array<Point> trianglePoints;
00681   Array<double> triangleWeights;
00682   getTriangleQuadPoints( trianglePoints , triangleWeights);
00683   double summWeights = 0.0;
00684 
00685   int nr1DPoints = linePoints.size();
00686   int nr2DPoints = quadQuadPoints.size();
00687 
00688   //those must be equal
00689   TEUCHOS_TEST_FOR_EXCEPTION( quadPoints.size() != quadQuadPoints.size() , std::runtime_error ,
00690       "quadPoints.size() != quadQuadPoints.size() , size1:" << quadPoints.size() << " , size2:" << quadQuadPoints.size());
00691   for (int q = 0; q < nr2DPoints; q++) {
00692     SUNDANCE_MSG3(verb_, " Quad point quadWeights["<<q<<"]="<<quadWeights[q]);
00693     summWeights = summWeights + quadWeights[q];
00694     quadWeights[q] = 0.0;
00695   }
00696   SUNDANCE_MSG3(verb_, " Summ old weights = " << summWeights );
00697 
00698     // for all elements in the list make the quadrature
00699   for (int listI = 0 ; listI < stackIndex ; listI++){
00700     // look at this quad only when it is not refined
00701     if (refinedStack[listI] < 0 ){
00702       SUNDANCE_MSG3(verb_, tabs << "getAdaptedQuadWeights Integrate quad listI:" << listI <<
00703           ",intersectioncaseStack[listI]:" << intersectioncaseStack[listI]);
00704       // ========== calculate first the integral over the whole quad
00705       Array<double> tmpWeightsQuad( quadWeights.size() , 0.0 );
00706       double  ofx , ofy , px , py;
00707       // integrate the whole quad, which result will be later used
00708             ofx = (quadPointStack[listI][1][0] - quadPointStack[listI][0][0]);
00709             ofy = (quadPointStack[listI][2][1] - quadPointStack[listI][0][1]);
00710             px = quadPointStack[listI][0][0]; py = quadPointStack[listI][0][1];
00711             // call the function which makes the quadrature
00712             makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00713                              quadQuadPoints , quadQuadWeights ,tmpWeightsQuad , 1.0 );
00714             SUNDANCE_MSG3(verb_, tabs << " end quadring the whole quad, now quadring the remaining parts " );
00715       switch (intersectioncaseStack[listI]){
00716       // =================================================================================================================
00717       case 0:{
00718         Array<double> tmpWeightsTriangle( quadWeights.size() , 0.0 );
00719         // integrate the triangle
00720                 ofx = (pointStack[listI][0][0] - quadPointStack[listI][0][0]);
00721                 ofy = (pointStack[listI][1][1] - quadPointStack[listI][0][1]);
00722                 px = quadPointStack[listI][0][0]; py = quadPointStack[listI][0][1];
00723                 // call the function which makes the quadrature
00724               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00725                   trianglePoints , triangleWeights ,tmpWeightsTriangle , 2.0 );
00726         // now add up the two diferent arreas
00727         for (int i = 0 ; i < nr2DPoints ; i++){
00728           SUNDANCE_MSG3(verb_, tabs << "i=" << i << ",Wquad[i]:" << tmpWeightsQuad[i] << " , Wtriag[i]:"<< tmpWeightsTriangle[i] );
00729           quadWeights[i] = alpha1[listI]*tmpWeightsTriangle[i] + alpha2[listI]*( tmpWeightsQuad[i] - tmpWeightsTriangle[i] );
00730         }
00731           break;}
00732       // =================================================================================================================
00733       case 1:{
00734         Array<double> tmpWeightsTriangle( quadWeights.size() , 0.0 );
00735         // integrate the triangle
00736                 ofx = (pointStack[listI][1][0] - quadPointStack[listI][0][0]);
00737                 ofy = (pointStack[listI][0][1] - quadPointStack[listI][2][1]);
00738                 px = quadPointStack[listI][2][0]; py = quadPointStack[listI][2][1];
00739                 // call the function which makes the quadrature
00740               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00741                   trianglePoints , triangleWeights ,tmpWeightsTriangle , 2.0 );
00742         // now add up the two different areas
00743         for (int i = 0 ; i < nr2DPoints ; i++){
00744           SUNDANCE_MSG3(verb_, tabs << "i=" << i << ",Wquad[i]:" << tmpWeightsQuad[i] << " , Wtriag[i]:"<< tmpWeightsTriangle[i] );
00745           quadWeights[i] = alpha2[listI]*tmpWeightsTriangle[i] + alpha1[listI]*( tmpWeightsQuad[i] - tmpWeightsTriangle[i] );
00746         }
00747           break;}
00748       // =================================================================================================================
00749       case 2:{
00750         Array<double> tmpWeightsTriangle( quadWeights.size() , 0.0 );
00751         // integrate the triangle
00752                 ofx = (pointStack[listI][0][0] - quadPointStack[listI][1][0]);
00753                 ofy = (pointStack[listI][1][1] - quadPointStack[listI][1][1]);
00754                 px = quadPointStack[listI][1][0]; py = quadPointStack[listI][1][1];
00755                 // call the function which makes the quadrature
00756               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00757                   trianglePoints , triangleWeights ,tmpWeightsTriangle , 2.0 );
00758         // now add up the two diferent arreas
00759         for (int i = 0 ; i < nr2DPoints ; i++){
00760           SUNDANCE_MSG3(verb_, tabs << "i=" << i << ",Wquad[i]:" << tmpWeightsQuad[i] << " , Wtriag[i]:"<< tmpWeightsTriangle[i] );
00761           quadWeights[i] = alpha2[listI]*tmpWeightsTriangle[i] + alpha1[listI]*( tmpWeightsQuad[i] - tmpWeightsTriangle[i] );
00762         }
00763           break;}
00764       // =================================================================================================================
00765       case 3:{
00766         Array<double> tmpWeightsTriangle( quadWeights.size() , 0.0 );
00767         // integrate the triangle
00768                 ofx = (pointStack[listI][1][0] - quadPointStack[listI][3][0]);
00769                 ofy = (pointStack[listI][0][1] - quadPointStack[listI][3][1]);
00770                 px = quadPointStack[listI][3][0]; py = quadPointStack[listI][3][1];
00771                 // call the function which makes the quadrature
00772               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00773                   trianglePoints , triangleWeights ,tmpWeightsTriangle , 2.0 );
00774         // now add up the two diferent arreas
00775         for (int i = 0 ; i < nr2DPoints ; i++){
00776           SUNDANCE_MSG3(verb_, tabs << "i=" << i << ",Wquad[i]:" << tmpWeightsQuad[i] << " , Wtriag[i]:"<< tmpWeightsTriangle[i] );
00777           quadWeights[i] = alpha2[listI]*tmpWeightsTriangle[i] + alpha1[listI]*( tmpWeightsQuad[i] - tmpWeightsTriangle[i] );
00778         }
00779           break;}
00780       // =================================================================================================================
00781       case 41: case 42:{
00782         // integrate the quad
00783         Array<double> tmpWeightsQuad2( quadWeights.size() , 0.0 );
00784         if (intersectioncaseStack[listI] == 41){
00785           ofx = ( pointStack[listI][1][0] - quadPointStack[listI][0][0]);
00786         } else {
00787           ofx = ( pointStack[listI][0][0] - quadPointStack[listI][0][0]);
00788         }
00789               ofy = (quadPointStack[listI][2][1] - quadPointStack[listI][0][1]);
00790               px = quadPointStack[listI][0][0]; py = quadPointStack[listI][0][1];
00791               // call the function which makes the quadrature
00792               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00793                   quadQuadPoints , quadQuadWeights ,tmpWeightsQuad2 , 1.0 );
00794         // integrate the triangle
00795         Array<double> tmpWeightsTriangle( quadWeights.size() , 0.0 );
00796         if (intersectioncaseStack[listI] == 41){
00797           ofx = ( pointStack[listI][0][0] - pointStack[listI][1][0]);
00798                   ofy = ( quadPointStack[listI][3][1] - quadPointStack[listI][0][1] );
00799           px = pointStack[listI][1][0]; py = pointStack[listI][0][1];
00800         } else {
00801           ofx = ( pointStack[listI][1][0] - pointStack[listI][0][0]);
00802           ofy = -( quadPointStack[listI][3][1] - quadPointStack[listI][0][1] );
00803           px = pointStack[listI][0][0]; py = pointStack[listI][1][1];
00804         }
00805         // call the function which makes the quadrature
00806               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00807                   trianglePoints , triangleWeights ,tmpWeightsTriangle , 2.0 );
00808         // now add up the two different areas
00809         for (int i = 0 ; i < nr2DPoints ; i++){
00810           SUNDANCE_MSG3(verb_, tabs << "i=" << i << " , Wquad[i]:" << tmpWeightsQuad[i] <<
00811               " , Wquad2[i]:" << tmpWeightsQuad2[i] << " , Wtriag[i]:"<< tmpWeightsTriangle[i] );
00812           quadWeights[i] = alpha1[listI]*( tmpWeightsQuad2[i] + tmpWeightsTriangle[i] ) +
00813                        alpha2[listI]*( tmpWeightsQuad[i] - tmpWeightsQuad2[i] - tmpWeightsTriangle[i]);
00814         }
00815           break;}
00816       // =================================================================================================================
00817       case 51: case 52:{
00818         // integrate the quad
00819         Array<double> tmpWeightsQuad2( quadWeights.size() , 0.0 );
00820         if (intersectioncaseStack[listI] == 52){
00821           ofy = ( pointStack[listI][0][1] - quadPointStack[listI][0][1]);
00822         } else {
00823           ofy = ( pointStack[listI][1][1] - quadPointStack[listI][0][1]);
00824         }
00825               ofx = (quadPointStack[listI][1][0] - quadPointStack[listI][0][0]);
00826               px = quadPointStack[listI][0][0]; py = quadPointStack[listI][0][1];
00827               // call the function which makes the quadrature
00828               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00829                   quadQuadPoints , quadQuadWeights ,tmpWeightsQuad2 , 1.0 );
00830         // integrate the triangle
00831         Array<double> tmpWeightsTriangle( quadWeights.size() , 0.0 );
00832         if (intersectioncaseStack[listI] == 51){
00833                   ofx = ( quadPointStack[listI][3][0] - quadPointStack[listI][0][0] );
00834           ofy = ( pointStack[listI][0][1] - pointStack[listI][1][1]);
00835           px = pointStack[listI][0][0]; py = pointStack[listI][1][1];
00836         } else {
00837           ofx = -( quadPointStack[listI][3][0] - quadPointStack[listI][0][0] );
00838           ofy = ( pointStack[listI][1][1] - pointStack[listI][0][1]);
00839           px = pointStack[listI][1][0]; py = pointStack[listI][0][1];
00840         }
00841         // call the function which makes the quadrature
00842               makeInterpolantQuad( px, py, ofx, ofy, nr1DPoints , nr2DPoints ,linePoints ,
00843                   trianglePoints , triangleWeights ,tmpWeightsTriangle , 2.0 );
00844         // now add up the two different areas
00845         for (int i = 0 ; i < nr2DPoints ; i++){
00846           SUNDANCE_MSG3(verb_, tabs << "i=" << i << " , Wquad[i]:" << tmpWeightsQuad[i] <<
00847               " , Wquad2[i]:" << tmpWeightsQuad2[i] << " , Wtriag[i]:"<< tmpWeightsTriangle[i] );
00848           quadWeights[i] = alpha1[listI]*(tmpWeightsQuad2[i]+tmpWeightsTriangle[i]) +
00849                        alpha2[listI]*(tmpWeightsQuad[i] - tmpWeightsQuad2[i] - tmpWeightsTriangle[i]);
00850         }
00851           break;}
00852       // =================================================================================================================
00853       case 6:{
00854         // no intersection point just integrate the quad elem
00855         for (int i = 0 ; i < nr2DPoints ; i++){
00856            SUNDANCE_MSG3(verb_, tabs << "i=" << i << " , Wquad[i]:" << tmpWeightsQuad[i] );
00857                    quadWeights[i] = quadWeights[i] + alpha1[listI]*tmpWeightsQuad[i];
00858         }
00859           break;}
00860       default:{
00861         // throw error
00862         TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error  , "Quad cell not integrable:" << intersectioncaseStack[listI]);
00863           break; }
00864       }
00865     }
00866   }
00867   // just print the weights
00868   summWeights = 0.0;
00869   for (int q = 0; q < nr2DPoints; q++) {
00870     summWeights = summWeights + quadWeights[q];
00871     SUNDANCE_MSG3(verb_, " New weights quadWeights["<<q<<"]="<<quadWeights[q]);
00872   }
00873   SUNDANCE_MSG3(verb_, " Summ new weights = " << summWeights );
00874 }
00875 
00876 
00877 // =========================================================================================================================
00878 
00879 void GaussLobattoQuadrature::getAdaptedQuadWeights_surf(
00880         int cellLID, const Mesh& mesh,
00881       const ParametrizedCurve& globalCurve, Array<Point>& quadPoints,
00882       Array<double>& quadWeights, bool& weightsChanged) const {
00883 
00884 
00885   int verb = 1;
00886   int nr_point = mesh.numFacets( mesh.spatialDim() , 0 , 0);
00887   Array<Point> maxCellsPoints(nr_point);
00888   int tmp_i , point_LID;
00889   weightsChanged = true;
00890 
00891   // - get all the brick cells point
00892   SUNDANCE_MSG3(verb, "CurveIntegralCalc::getCurveQuadPoints nr points per cell: " << nr_point)
00893   for (int jj = 0 ; jj < nr_point ; jj++){
00894     point_LID = mesh.facetLID( mesh.spatialDim() , cellLID , 0 , jj , tmp_i );
00895     maxCellsPoints[jj] = mesh.nodePosition(point_LID);
00896     SUNDANCE_MSG3(verb, " max cell point p[" << jj << "]:"<< maxCellsPoints[jj]);
00897   }
00898 
00899   // - use the existing method to get the triangles
00900   Array<Point> intersectPoints;
00901   Array<int> edgeIndex;
00902   Array<int> triangleIndex;
00903     // the results in this form will be in physical coordinates
00904   SundanceSurf3DCalc::getCurveQuadPoints( BrickCell , cellLID , mesh , globalCurve, maxCellsPoints,
00905              intersectPoints ,  edgeIndex , triangleIndex );
00906 
00907   // initialize the output variables
00908   getBrickRule(quadPoints,quadWeights);
00909 
00910   // if the cell has no intersection points then just return the standard weights
00911   if (intersectPoints.size() < 1){
00912     double intCoeff = globalCurve.integrationParameter( 0.5*(maxCellsPoints[0]+maxCellsPoints[7]) );
00913     SUNDANCE_MSG1(verb, " WARNING cel has not been intersected , alpha = " << intCoeff );
00914     for (int ii = 0 ; ii < quadWeights.size() ; ii++ ) { quadWeights[ii] = quadWeights[ii] * intCoeff;}
00915     return;
00916   }
00917 
00918   // - determine the projection direction , and the points which overlap with the projected quad cell
00919   int nrIntPoint = intersectPoints.size() , projDir = -1 , tmpI , coveredNodes = 0 , firstCoveredNode = -1;
00920   int possibleProject[3] = { -1 , -1 , -1};
00921   int otherDims[2] = {-1,-1};
00922   int projectedEdgeCovered[4] = { -1 , -1 , -1 , -1 };
00923   double dist_tmp = 1e+100;
00924 
00925   // find out the possible projections (there is only three possible projections)
00926   for (int ii = 0 ; ii < nrIntPoint ; ii++){
00927     // for each edge index we set the projection directory to one
00928     possibleProject[ edge3DProjection[ edgeIndex[ii] ] ] = 1;
00929   }
00930   // here we find the projection directory
00931   for (int ii = 0 ; ii < 3 ; ii++){
00932     if (possibleProject[ii] > 0) {
00933       SUNDANCE_MSG3(verb, " Projected dir = " << ii);
00934       projDir = ii; break;
00935     }
00936   }
00937   // here we store the two other directions which are not projected
00938   tmpI = 0;
00939   for (int ii = 0 ; ii < 3 ; ii++){
00940     if ( ii != projDir ) {
00941       SUNDANCE_MSG3(verb, " otherDims["<< tmpI << "] = " << ii);
00942       otherDims[tmpI] = ii; tmpI++;
00943     }
00944   }
00945 
00946   // - transform all intersection points to reference coordinates
00947   for (int ii = 0 ; ii < nrIntPoint ; ii++){
00948     intersectPoints[ii] = Point( (intersectPoints[ii][0] - maxCellsPoints[0][0])/(maxCellsPoints[7][0] - maxCellsPoints[0][0]) ,
00949          (intersectPoints[ii][1] - maxCellsPoints[0][1])/(maxCellsPoints[7][1] - maxCellsPoints[0][1]) ,
00950          (intersectPoints[ii][2] - maxCellsPoints[0][2])/(maxCellsPoints[7][2] - maxCellsPoints[0][2]));
00951     SUNDANCE_MSG3(verb, "REF intersectPoints[" << ii << "] = "<< intersectPoints[ii]);
00952   }
00953 
00954   // - the quad points on the projected
00955   Array<Point> quadProjectPoints(4);
00956   quadProjectPoints[0] = Point(0.0,0.0); quadProjectPoints[1] = Point(1.0,0.0);
00957   quadProjectPoints[2] = Point(0.0,1.0); quadProjectPoints[3] = Point(1.0,1.0);
00958 
00959   // - determine which out of the 4 nodes has been covered
00960   coveredNodes = 0;
00961   for (int ii = 0 ; ii < nrIntPoint ; ii++) {
00962     Point tmpP( intersectPoints[ii][otherDims[0]] , intersectPoints[ii][otherDims[1]] );
00963     for (int tt = 0 ; tt < 4 ; tt++ ) {
00964       if ( ((tmpP-quadProjectPoints[tt])*(tmpP-quadProjectPoints[tt]) < 1e-12) && (projectedEdgeCovered[tt] < 0) )
00965       {
00966         projectedEdgeCovered[tt] = 1; // this node is covered
00967         //SUNDANCE_MSG3(verb, "COVERED node[" << tt << "] = "<< quadProjectPoints[tt]);
00968         coveredNodes++;
00969         // to store the first node which is covered is important to measure the integration coefficient
00970         if (firstCoveredNode < 0) { firstCoveredNode = tt; }
00971       }
00972     }
00973   }
00974   SUNDANCE_MSG3( verb , "covered nodes = "<< coveredNodes );
00975 
00976   // if not all nodes in the projected quad are covered then we have to add them
00977   // so we loop as long this is done
00978   while (coveredNodes < 4) {
00979       // - look for one edge, which has at least one neighbor which is covered
00980     int foundEdgeIndex[3] = {-1,-1,-1} , foundI = 0 , pF = -1 , pL = -1;
00981     int neighbourNodes[4][2] = { {1,2} , {0,3} , {0,3} , {1,2} };
00982     for (int ii = 0 ; ii < 4 ; ii++){
00983       //SUNDANCE_MSG3( verb , "test node = "<< ii << " , projectedEdgeCovered[ii]=" << projectedEdgeCovered[ii] <<
00984       //    " , N1=" << projectedEdgeCovered[ neighbourNodes[ii][0] ] << " , N2=" << projectedEdgeCovered[ neighbourNodes[ii][1] ]);
00985       // look for uncovered nodes, which have one neighbor which is covered
00986       if ( (projectedEdgeCovered[ii] < 0) &&
00987          ((projectedEdgeCovered[ neighbourNodes[ii][0] ] > 0) || (projectedEdgeCovered[ neighbourNodes[ii][1] ] > 0)) ){
00988         //SUNDANCE_MSG3( verb , "found node = "<< ii );
00989         foundEdgeIndex[foundI] = ii;
00990         projectedEdgeCovered[ii] = 1;
00991         coveredNodes++;
00992         break;
00993       }
00994     }
00995     foundI = 0;
00996     // look for the nearest point, with respect to the found node
00997     dist_tmp = 1e+100;
00998     for (int ii = 0 ; ii < nrIntPoint ; ii++){
00999       Point tmpP( intersectPoints[ii][otherDims[0]] , intersectPoints[ii][otherDims[1]] );
01000       double dist_tmp_tmp = ::sqrt( (tmpP-quadProjectPoints[foundEdgeIndex[foundI]])*(tmpP-quadProjectPoints[foundEdgeIndex[foundI]]) );
01001       //SUNDANCE_MSG3( verb , "test projected Int Point = "<< tmpP << " , dist = " << dist_tmp_tmp);
01002       if ( dist_tmp_tmp < dist_tmp) { pF = ii; dist_tmp = dist_tmp_tmp;}
01003     }
01004     SUNDANCE_MSG3( verb , "pF = "<< pF );
01005 
01006     // loop as long all the neighbors of the neighbors are covered
01007     while ( (projectedEdgeCovered[ neighbourNodes[foundEdgeIndex[foundI]][0] ] < 0) ||
01008           (projectedEdgeCovered[ neighbourNodes[foundEdgeIndex[foundI]][1] ] < 0) ) {
01009       int tmp_next_Neighbor = -1;
01010       // only one can be true
01011       if (projectedEdgeCovered[ neighbourNodes[foundEdgeIndex[foundI]][0] ] < 0) { tmp_next_Neighbor = neighbourNodes[foundEdgeIndex[foundI]][0]; }
01012       if (projectedEdgeCovered[ neighbourNodes[foundEdgeIndex[foundI]][1] ] < 0) { tmp_next_Neighbor = neighbourNodes[foundEdgeIndex[foundI]][1]; }
01013       //SUNDANCE_MSG3( verb , "tmp_next_Neighbor = "<< tmp_next_Neighbor );
01014       foundEdgeIndex[foundI+1] = tmp_next_Neighbor;
01015       projectedEdgeCovered[tmp_next_Neighbor] = 1;
01016       foundI++;
01017       coveredNodes++;
01018     }
01019     // once we finish the loop we look for a next intersection point which is closer to the last node and is different from "pF"
01020     dist_tmp = 1e+100;
01021     for (int ii = 0 ; ii < nrIntPoint ; ii++){
01022       Point tmpP( intersectPoints[ii][otherDims[0]] , intersectPoints[ii][otherDims[1]] );
01023       double dist_tmp_tmp = ::sqrt((tmpP-quadProjectPoints[foundEdgeIndex[foundI]])*(tmpP-quadProjectPoints[foundEdgeIndex[foundI]]));
01024       if ( (dist_tmp_tmp < dist_tmp) && (ii != pF ) ) { pL = ii; dist_tmp = dist_tmp_tmp; }
01025     }
01026     SUNDANCE_MSG3( verb , "pL = "<< pF );
01027 
01028     // - now we have a list of the nodes of the projected quad, which need to be added
01029     for (int nn = 0 ; nn <= foundI ; nn = nn + 1)
01030     {
01031       int firstI_Tri = pF, secondI_Tri = -1 , thirdI_Tri = -1;
01032       intersectPoints.resize( nrIntPoint + 1);
01033       nrIntPoint++;
01034       //  check if one point will not be added twice ... on the other side this is not a problem for the integration
01035       Point ptmp( 0 , 0 , 0 );
01036       ptmp[otherDims[0]] = quadProjectPoints[foundEdgeIndex[nn]][0];
01037       ptmp[otherDims[1]] = quadProjectPoints[foundEdgeIndex[nn]][1];
01038       ptmp[projDir] = intersectPoints[pF][projDir];
01039       //SUNDANCE_MSG3( verb , "Add Point 2D " << quadProjectPoints[foundEdgeIndex[nn]] << " , P:" << ptmp);
01040       intersectPoints[nrIntPoint-1] = ptmp;
01041       secondI_Tri = nrIntPoint -1;
01042       Point nextP( 0 , 0 , 0 );
01043       if (nn < foundI ){
01044         nextP[otherDims[0]] = quadProjectPoints[foundEdgeIndex[nn+1]][0];
01045         nextP[otherDims[1]] = quadProjectPoints[foundEdgeIndex[nn+1]][1];
01046         nextP[projDir] = intersectPoints[pF][projDir];
01047         //SUNDANCE_MSG3( verb , "Add Point 2D " << quadProjectPoints[foundEdgeIndex[nn+1]] );
01048         intersectPoints.resize( nrIntPoint + 1);
01049         nrIntPoint++;
01050         thirdI_Tri = nrIntPoint - 1;
01051         intersectPoints[nrIntPoint-1] = nextP;
01052       } else {
01053         // the last point is "pL"
01054         nextP = intersectPoints[pL];
01055         thirdI_Tri = pL;
01056       }
01057       int triagSize = triangleIndex.size();
01058       triangleIndex.resize(triagSize + 3);
01059       triangleIndex[triagSize] = firstI_Tri;
01060       triangleIndex[triagSize+1] = secondI_Tri;
01061       triangleIndex[triagSize+2] = thirdI_Tri;
01062       SUNDANCE_MSG3( verb , "Add triangle [ "<< firstI_Tri << " , " <<  secondI_Tri << " , " << thirdI_Tri << " ]");
01063       SUNDANCE_MSG3( verb , "Added triangle's points [ "<< intersectPoints[firstI_Tri] << " , "
01064           <<  intersectPoints[secondI_Tri] << " , " << intersectPoints[thirdI_Tri] << " ]");
01065     }
01066 
01067     SUNDANCE_MSG3( verb , "covered nodes = "<< coveredNodes );
01068   } //end from the while which makes sure that each nodes are covered
01069 
01070 
01071   // - first get the quadrature points for the whole quad
01072   Array<Point> wholeQuadPoints;
01073   Array<double> wholeQuadweights;
01074   getBrickRule(wholeQuadPoints,wholeQuadweights);
01075   Array<Point> linePoints;
01076   Array<double> lineWeights;
01077   getLineRule(linePoints,lineWeights);
01078   int nrLinePoints = linePoints.size() , nr3DPoints = wholeQuadPoints.size();
01079   // this will be the output of the quadrature
01080   Array<double> computedWeights(nr3DPoints,0.0);
01081 
01082   // compute the quadrature points and weights for the prism elements
01083   Array<Point> triagPoints;
01084   Array<double> triagWeights;
01085   getTriangleQuadPoints( triagPoints , triagWeights );
01086   int nrQuadTriagPoints = triagPoints.size();
01087   int nrQuadPrism = triagPoints.size() * linePoints.size();
01088   Array<Point> prismPoints(nrQuadPrism);
01089   Array<double> prismWeights(nrQuadPrism);
01090 
01091   // - integrate all the prism elements, if the heights are different then zero
01092   int nr_prism = triangleIndex.size()/3;
01093   for (int p = 0 ; p < nr_prism ; p++ ){
01094     Point p0 = intersectPoints[ triangleIndex[3*p] ];
01095     Point p1 = intersectPoints[ triangleIndex[3*p + 1] ];
01096     Point p2 = intersectPoints[ triangleIndex[3*p + 2] ];
01097     Point p0L = p0;
01098     Point p1L = p1;
01099     Point p2L = p2;
01100     p0L[projDir] = 0.0; p1L[projDir] = 0.0; p2L[projDir] = 0.0;
01101     //SUNDANCE_MSG3( verb , " Triangle [ " << p0 << " , "  << p1 << " , " << p2 << "]");
01102     //SUNDANCE_MSG3( verb , " Triangle bottom [ " << p0L << " , "  << p1L << " , " << p2L << "]");
01103 
01104     // only integrate the prism when there is a measurable height
01105     if ( (::fabs(p0[projDir]-p0L[projDir]) + ::fabs(p1[projDir]-p1L[projDir])+ ::fabs(p2[projDir]-p2L[projDir]) ) > 1e-12 ){
01106       Point normAreaVect = cross(p1L-p0L ,p2L-p0L);
01107       Point triagPointTmp;
01108       double areaTriag = 0.5*::sqrt(normAreaVect*normAreaVect) , Zc = 0.0;
01109       int quadPIndex = 0;
01110       // compute the quadrature points for this prism
01111       for (int triagQIndex = 0 ; triagQIndex < nrQuadTriagPoints ; triagQIndex++)
01112       {
01113         // Zc is the height of the prism at the quadrature position
01114           Zc = (1.0 - triagPoints[triagQIndex][0] - triagPoints[triagQIndex][1]) * ( p0[projDir] - p0L[projDir] );
01115         Zc = Zc + triagPoints[triagQIndex][0]*( p1[projDir] - p1L[projDir] );
01116         Zc = Zc + triagPoints[triagQIndex][1]*( p2[projDir] - p2L[projDir] );
01117         triagPointTmp = p0L + triagPoints[triagQIndex][0]*(p1L-p0L) + triagPoints[triagQIndex][1]*(p2L-p0L);
01118         //SUNDANCE_MSG3( verb , " areaTriag = " << areaTriag << " , Zc = "  << Zc <<
01119         //    " , triagWeights[triagQIndex] = " << triagWeights[triagQIndex] << " , lineWeights[0]" << lineWeights[0]);
01120         for (int lineQIndex = 0 ; lineQIndex < nrLinePoints ; lineQIndex++ , quadPIndex++)
01121         {
01122           // take the tensor product of the triangle and the line quadrature points
01123           // it is important that we compute correctly the weights and the quadrature points positions
01124           prismPoints[quadPIndex] = triagPointTmp;
01125           prismPoints[quadPIndex][projDir] = Zc * linePoints[lineQIndex][0];
01126           prismWeights[quadPIndex] = Zc*areaTriag*triagWeights[triagQIndex]*lineWeights[lineQIndex];
01127         }
01128       }
01129       // call the method to integrate the Lagrange polynoms
01130       makeInterpolantBrick( nrLinePoints , nr3DPoints , linePoints ,
01131           prismPoints , prismWeights , computedWeights , 1.0);
01132     }
01133   }
01134 
01135   // - determine the integration coefficient , use "firstCoveredNode" with "0" and "1" offset
01136   Point tmpPointIntCoeff = intersectPoints[firstCoveredNode];
01137   tmpPointIntCoeff[projDir] = -0.001;
01138   Point tmpPointIntCoeff1 = Point( maxCellsPoints[0][0] + tmpPointIntCoeff[0]*(maxCellsPoints[7][0]-maxCellsPoints[0][0]) ,
01139                                maxCellsPoints[0][1] + tmpPointIntCoeff[1]*(maxCellsPoints[7][1]-maxCellsPoints[0][1]) ,
01140                                maxCellsPoints[0][2] + tmpPointIntCoeff[2]*(maxCellsPoints[7][2]-maxCellsPoints[0][2]) );
01141   double zeroIntCoeff = globalCurve.integrationParameter( tmpPointIntCoeff1 );
01142   tmpPointIntCoeff[projDir] = 1.001;
01143   tmpPointIntCoeff1 = Point( maxCellsPoints[0][0] + tmpPointIntCoeff[0]*(maxCellsPoints[7][0]-maxCellsPoints[0][0]) ,
01144                        maxCellsPoints[0][1] + tmpPointIntCoeff[1]*(maxCellsPoints[7][1]-maxCellsPoints[0][1]) ,
01145                        maxCellsPoints[0][2] + tmpPointIntCoeff[2]*(maxCellsPoints[7][2]-maxCellsPoints[0][2]) );
01146   double oneIntCoeff = globalCurve.integrationParameter( tmpPointIntCoeff1 );
01147 
01148   SUNDANCE_MSG3( verb , " alpha1 = " << zeroIntCoeff << " , alpha2 = "  << oneIntCoeff );
01149 
01150   // at the end just compute the weights, by simply taking the difference of the whole and the computed weights
01151   double sumWeights = 0.0;
01152   for (int q = 0 ; q < computedWeights.size() ; q++ ){
01153     quadWeights[q] = zeroIntCoeff * computedWeights[q] + oneIntCoeff*( wholeQuadweights[q]- computedWeights[q]);
01154     sumWeights = sumWeights + quadWeights[q];
01155   }
01156 
01157   SUNDANCE_MSG3(verb , " ---------- END METHOD -----------  sum weights = " << sumWeights);
01158 
01159 }
01160 
01161 
01162 // ======================================================================================================
01163 
01164 
01165 
01166 
01167 void GaussLobattoQuadrature::getLineRule(
01168     Array<Point>& quadPointsL,
01169     Array<double>& quadWeights) const {
01170 
01171   int n = order()+1; //this is correct m points for a basis of order m-1
01172 
01173   Array<double> quadPoints;
01174   quadPoints.resize(n);
01175   quadPointsL.resize(n);
01176   quadWeights.resize(n);
01177 
01178   // ================== QUAD POINTS ========================
01179   switch (n) {
01180   case 2: { quadPoints[0]=0.0; quadPoints[1] = 1; break; }
01181   case 3: { quadPoints[0]=0.0; quadPoints[1] = 0.5; quadPoints[2] = 1.0; break; }
01182   case 4: { quadPoints[0]=0.0; quadPoints[1] = 0.5-0.5/::sqrt(5.0); quadPoints[2] = 0.5+0.5/::sqrt(5.0); quadPoints[3] = 1.0; break; }
01183   case 5: { quadPoints[0]=0.0; quadPoints[1] = 0.5-0.5*::sqrt(3.0/7.0); quadPoints[2] = 0.5;
01184           quadPoints[3] = 0.5+0.5*::sqrt(3.0/7.0); quadPoints[4] = 1.0; break; }
01185   case 6: { double t0=::sqrt(7.0);
01186           double t1=(7.0+2.0*t0)/21.0;
01187           double t2=(7.0-2.0*t0)/21.0;
01188           quadPoints[0] = 0; quadPoints[1] = 0.5-0.5*::sqrt(t1); quadPoints[2] = 0.5-0.5*::sqrt(t2);
01189           quadPoints[3] = 0.5+0.5*::sqrt(t2); quadPoints[4] = 0.5+0.5*::sqrt(t1); quadPoints[5] = 1.0; break; }
01190   case 7: {
01191     quadPoints[0] = 0.00000000000000000000e+00;
01192     quadPoints[1] = 8.48880518607165179823e-02;
01193     quadPoints[2] = 2.65575603264642912116e-01;
01194     quadPoints[3] = 5.00000000000000000000e-01;
01195     quadPoints[4] = 7.34424396735357087884e-01;
01196     quadPoints[5] = 9.15111948139283537529e-01;
01197     quadPoints[6] = 1.00000000000000000000e+00;
01198       break; }
01199   case 8: {
01200     quadPoints[0] = 0.00000000000000000000e+00;
01201     quadPoints[1] = 6.41299257451967141819e-02;
01202     quadPoints[2] = 2.04149909283428854234e-01;
01203     quadPoints[3] = 3.95350391048760574364e-01;
01204     quadPoints[4] = 6.04649608951239425636e-01;
01205     quadPoints[5] = 7.95850090716571090255e-01;
01206     quadPoints[6] = 9.35870074254803285818e-01;
01207     quadPoints[7] = 1.00000000000000000000e+00;
01208     break; }
01209   case 9: {
01210     quadPoints[0] = 0.00000000000000000000e+00;
01211     quadPoints[1] = 5.01210022942699118254e-02;
01212     quadPoints[2] = 1.61406860244631134016e-01;
01213     quadPoints[3] = 3.18441268086910922452e-01;
01214     quadPoints[4] = 5.00000000000000000000e-01;
01215     quadPoints[5] = 6.81558731913089133059e-01;
01216     quadPoints[6] = 8.38593139755368865984e-01;
01217     quadPoints[7] = 9.49878997705730032663e-01;
01218     quadPoints[8] = 1.00000000000000000000e+00;
01219       break; }
01220   case 10: {
01221     quadPoints[0] = 0.00000000000000000000e+00;
01222     quadPoints[1] = 4.02330459167705711820e-02;
01223     quadPoints[2] = 1.30613067447247432895e-01;
01224     quadPoints[3] = 2.61037525094777733692e-01;
01225     quadPoints[4] = 4.17360521166806497373e-01;
01226     quadPoints[5] = 5.82639478833193447116e-01;
01227     quadPoints[6] = 7.38962474905222266308e-01;
01228     quadPoints[7] = 8.69386932552752567105e-01;
01229     quadPoints[8] = 9.59766954083229428818e-01;
01230     quadPoints[9] = 1.00000000000000000000e+00;
01231     break; }
01232    case 11: {
01233     quadPoints[0] = 0.00000000000000000000e+00;
01234     quadPoints[1] = 3.29992847959704183047e-02;
01235     quadPoints[2] = 1.07758263168427792511e-01;
01236     quadPoints[3] = 2.17382336501897477365e-01;
01237     quadPoints[4] = 3.52120932206530290465e-01;
01238     quadPoints[5] = 5.00000000000000000000e-01;
01239     quadPoints[6] = 6.47879067793469709535e-01;
01240     quadPoints[7] = 7.82617663498102578146e-01;
01241     quadPoints[8] = 8.92241736831572263000e-01;
01242     quadPoints[9] = 9.67000715204029637206e-01;
01243     quadPoints[10] = 1.00000000000000000000e+00;
01244       break; }
01245    case 12: {
01246     quadPoints[0] = 0.00000000000000000000e+00;
01247     quadPoints[1] = 2.75503638885589152707e-02;
01248     quadPoints[2] = 9.03603391779966846897e-02;
01249     quadPoints[3] = 1.83561923484069688950e-01;
01250     quadPoints[4] = 3.00234529517325543502e-01;
01251     quadPoints[5] = 4.31723533572536233294e-01;
01252     quadPoints[6] = 5.68276466427463766706e-01;
01253     quadPoints[7] = 6.99765470482674456498e-01;
01254     quadPoints[8] = 8.16438076515930255539e-01;
01255     quadPoints[9] = 9.09639660822003315310e-01;
01256     quadPoints[10] = 9.72449636111441084729e-01;
01257     quadPoints[11] = 1.00000000000000000000e+00;
01258       break; }
01259   }
01260 
01261   // transform the array of doubles into an array of Points
01262   for (int i = 0 ; i < n ; i++){
01263     quadPointsL[i] = Point(quadPoints[i]);
01264   }
01265 
01266   // ================== WEIGHTS ========================
01267   switch (n) {
01268   case 2:{
01269     quadWeights[0] = 0.5;
01270     quadWeights[1] = 0.5;
01271       break; }
01272   case 3:{
01273     quadWeights[0] = 1.0/6.0; quadWeights[1] = 2.0/3.0; quadWeights[2] = 1.0/6.0;
01274     break;}
01275   case 4:{
01276     quadWeights[0] = 1.0/12.0; quadWeights[1] = 5.0/12.0; quadWeights[2] = 5.0/12.0; quadWeights[3] = 1.0/12.0;
01277     break;}
01278   case 5:{
01279     quadWeights[0] = 0.05; quadWeights[1] = 49.0/180.0; quadWeights[2] = 32.0/90.0; quadWeights[3] = 49.0/180.0; quadWeights[4] = 0.05;
01280     break;}
01281   case 6:{
01282       double t0=::sqrt(7.0);
01283       double t1=(7.0+2.0*t0)/21.0;
01284       double t2=(7.0-2.0*t0)/21.0;
01285       double k1=(1.0-t0)*(1.0-t0);
01286       double k2=(1.0+t0)*(1.0+t0);
01287       quadWeights[0] = 1.0/30.0; quadWeights[1] = 0.3/(t1*k1); quadWeights[2] = 0.3/(t2*k2); quadWeights[3] = 0.3/(t2*k2);
01288       quadWeights[4] = 0.3/(t1*k1); quadWeights[5] = 1.0/30.0;
01289       break;}
01290   case 7:{
01291     quadWeights[0] = 2.38095238095238082021e-02;
01292     quadWeights[1] = 1.38413023680782953928e-01;
01293     quadWeights[2] = 2.15872690604931305458e-01;
01294     quadWeights[3] = 2.43809523809523809312e-01;
01295     quadWeights[4] = 2.15872690604931305458e-01;
01296     quadWeights[5] = 1.38413023680782953928e-01;
01297     quadWeights[6] = 2.38095238095238082021e-02;
01298       break;}
01299   case 8:{
01300     quadWeights[0] = 1.78571428571428561516e-02;
01301     quadWeights[1] = 1.05352113571753072674e-01;
01302     quadWeights[2] = 1.70561346241752204156e-01;
01303     quadWeights[3] = 2.06229397329351860080e-01;
01304     quadWeights[4] = 2.06229397329351860080e-01;
01305     quadWeights[5] = 1.70561346241752204156e-01;
01306     quadWeights[6] = 1.05352113571753072674e-01;
01307     quadWeights[7] = 1.78571428571428561516e-02;
01308       break;}
01309   case 9:{
01310     quadWeights[0] = 1.38888888888888881179e-02;
01311     quadWeights[1] = 8.27476807804027880699e-02;
01312     quadWeights[2] = 1.37269356250080826198e-01;
01313     quadWeights[3] = 1.73214255486523083238e-01;
01314     quadWeights[4] = 1.85759637188208620584e-01;
01315     quadWeights[5] = 1.73214255486523083238e-01;
01316     quadWeights[6] = 1.37269356250080826198e-01;
01317     quadWeights[7] = 8.27476807804027880699e-02;
01318     quadWeights[8] = 1.38888888888888881179e-02;
01319       break;}
01320   case 10:{
01321     quadWeights[0] = 1.11111111111111115352e-02;
01322     quadWeights[1] = 6.66529954255350304271e-02;
01323     quadWeights[2] = 1.12444671031563220298e-01;
01324     quadWeights[3] = 1.46021341839841889421e-01;
01325     quadWeights[4] = 1.63769880591948718829e-01;
01326     quadWeights[5] = 1.63769880591948718829e-01;
01327     quadWeights[6] = 1.46021341839841889421e-01;
01328     quadWeights[7] = 1.12444671031563220298e-01;
01329     quadWeights[8] = 6.66529954255350304271e-02;
01330     quadWeights[9] = 1.11111111111111115352e-02;
01331       break;}
01332   case 11:{
01333     quadWeights[0] = 9.09090909090909046752e-03;
01334     quadWeights[1] = 5.48061366334974126024e-02;
01335     quadWeights[2] = 9.35849408901525958715e-02;
01336     quadWeights[3] = 1.24024052132014145355e-01;
01337     quadWeights[4] = 1.43439562389503921791e-01;
01338     quadWeights[5] = 1.50108797727845355574e-01;
01339     quadWeights[6] = 1.43439562389503921791e-01;
01340     quadWeights[7] = 1.24024052132014145355e-01;
01341     quadWeights[8] = 9.35849408901525958715e-02;
01342     quadWeights[9] = 5.48061366334974126024e-02;
01343     quadWeights[10] = 9.09090909090909046752e-03;
01344       break;}
01345   case 12:{
01346     quadWeights[0] = 7.57575757575757596785e-03;
01347     quadWeights[1] = 4.58422587065981240739e-02;
01348     quadWeights[2] = 7.89873527821850218711e-02;
01349     quadWeights[3] = 1.06254208880510653268e-01;
01350     quadWeights[4] = 1.25637801599600640312e-01;
01351     quadWeights[5] = 1.35702620455348088591e-01;
01352     quadWeights[6] = 1.35702620455348088591e-01;
01353     quadWeights[7] = 1.25637801599600640312e-01;
01354     quadWeights[8] = 1.06254208880510653268e-01;
01355     quadWeights[9] = 7.89873527821850218711e-02;
01356     quadWeights[10] = 4.58422587065981240739e-02;
01357     quadWeights[11] = 7.57575757575757596785e-03;
01358       break;}
01359    }
01360 }
01361 
01362 
01363 void GaussLobattoQuadrature::getTriangleQuadPoints(Array<Point>& pnt  ,Array<double>& weight ) const{
01364   // we took this directly from the Matlab code of Prof.Ulbrich
01365   int deg = nrPointin1D_ + nrPointin1D_ - 2;
01366   if (deg==2){
01367     pnt.resize(3); weight.resize(3);
01368     pnt[0] = Point(0.50000000000000000000 , 0.00000000000000000000);
01369     pnt[1] = Point(0.50000000000000000000 , 0.50000000000000000000);
01370     pnt[2] = Point(0.00000000000000000000 , 0.50000000000000000000);
01371     weight[0] = 0.33333333333333333333;
01372     weight[1] = 0.33333333333333333333;
01373     weight[2] = 0.33333333333333333333;
01374   } else if (deg==3){
01375     pnt.resize(4); weight.resize(4);
01376     pnt[0] = Point(0.33333333333333333333 , 0.33333333333333333333);
01377     pnt[1] = Point(0.60000000000000000000 , 0.20000000000000000000);
01378     pnt[2] = Point(0.20000000000000000000 , 0.60000000000000000000);
01379     pnt[3] = Point(0.20000000000000000000 , 0.20000000000000000000);
01380     weight[0] = -0.56250000000000000000;
01381     weight[1] = 0.52083333333333333333;
01382     weight[2] = 0.52083333333333333333;
01383     weight[3] = 0.52083333333333333333;
01384   } else if (deg==4) {
01385     pnt.resize(6); weight.resize(6);
01386     pnt[0] = Point(0.816847572980459 , 0.091576213509771);
01387     pnt[1] = Point(0.091576213509771 , 0.816847572980459);
01388     pnt[2] = Point(0.091576213509771 , 0.091576213509771);
01389     pnt[3] = Point(0.108103018168070 , 0.445948490915965);
01390     pnt[4] = Point(0.445948490915965 , 0.108103018168070);
01391     pnt[5] = Point(0.445948490915965 , 0.445948490915965);
01392     weight[0] = 0.109951743655322;
01393     weight[1] = 0.109951743655322;
01394     weight[2] = 0.109951743655322;
01395     weight[3] = 0.223381589678011;
01396     weight[4] = 0.223381589678011;
01397     weight[5] = 0.223381589678011;
01398   } else if (deg==5) {
01399     pnt.resize(7); weight.resize(7);
01400     pnt[0] = Point(0.33333333333333333 , 0.33333333333333333);
01401     pnt[1] = Point(0.79742698535308720 , 0.10128650732345633);
01402     pnt[2] = Point(0.10128650732345633 , 0.79742698535308720);
01403     pnt[3] = Point(0.10128650732345633 , 0.10128650732345633);
01404     pnt[4] = Point(0.05971587178976981 , 0.47014206410511505);
01405     pnt[5] = Point(0.47014206410511505 , 0.05971587178976981);
01406     pnt[6] = Point(0.47014206410511505 , 0.47014206410511505);
01407     weight[0] = 0.22500000000000000;
01408     weight[1] = 0.12593918054482717;
01409     weight[2] = 0.12593918054482717;
01410     weight[3] = 0.12593918054482717;
01411     weight[4] = 0.13239415278850616;
01412     weight[5] = 0.13239415278850616;
01413     weight[6] = 0.13239415278850616;
01414   } else if (deg==6) {
01415     pnt.resize(9); weight.resize(9);
01416     pnt[0] = Point(0.124949503233232 , 0.437525248383384);
01417     pnt[1] = Point(0.437525248383384 , 0.124949503233232);
01418     pnt[2] = Point(0.437525248383384 , 0.437525248383384);
01419     pnt[3] = Point(0.797112651860071 , 0.165409927389841);
01420     pnt[4] = Point(0.797112651860071 , 0.037477420750088);
01421     pnt[5] = Point(0.165409927389841 , 0.797112651860071);
01422     pnt[6] = Point(0.165409927389841 , 0.037477420750088);
01423     pnt[7] = Point(0.037477420750088 , 0.797112651860071);
01424     pnt[8] = Point(0.037477420750088 , 0.165409927389841);
01425     weight[0] = 0.205950504760887;
01426     weight[1] = 0.205950504760887;
01427     weight[2] = 0.205950504760887;
01428     weight[3] = 0.063691414286223;
01429     weight[4] = 0.063691414286223;
01430     weight[5] = 0.063691414286223;
01431     weight[6] = 0.063691414286223;
01432     weight[7] = 0.063691414286223;
01433     weight[8] = 0.063691414286223;
01434   } else if (deg==7) {
01435     pnt.resize(13); weight.resize(13);
01436     pnt[0] = Point(0.333333333333333 , 0.333333333333333);
01437     pnt[1] = Point(0.479308067841923 , 0.260345966079038);
01438     pnt[2] = Point(0.260345966079038 , 0.479308067841923);
01439     pnt[3] = Point(0.260345966079038 , 0.260345966079038);
01440     pnt[4] = Point(0.869739794195568 , 0.065130102902216);
01441     pnt[5] = Point(0.065130102902216 , 0.869739794195568);
01442     pnt[6] = Point(0.065130102902216  ,0.065130102902216);
01443     pnt[7] = Point(0.638444188569809 , 0.312865496004875);
01444     pnt[8] = Point(0.638444188569809 , 0.048690315425316);
01445     pnt[9] = Point(0.312865496004875 , 0.638444188569809);
01446     pnt[10] = Point(0.312865496004875 , 0.048690315425316);
01447     pnt[11] = Point(0.048690315425316 , 0.638444188569809);
01448     pnt[12] = Point(0.048690315425316 , 0.312865496004875);
01449     weight[0] = -0.149570044467670;
01450     weight[1] = 0.175615257433204;
01451     weight[2] = 0.175615257433204;
01452     weight[3] = 0.175615257433204;
01453     weight[4] = 0.053347235608839;
01454     weight[5] = 0.053347235608839;
01455     weight[6] = 0.053347235608839;
01456     weight[7] = 0.077113760890257;
01457     weight[8] = 0.077113760890257;
01458     weight[9] = 0.077113760890257;
01459     weight[10] = 0.077113760890257;
01460     weight[11] = 0.077113760890257;
01461     weight[12] = 0.077113760890257;
01462     } else if (deg==8) {
01463     pnt.resize(19); weight.resize(19);
01464     pnt[0] = Point(0.3333333333333333 , 0.3333333333333333);
01465     pnt[1] = Point(0.7974269853530872 , 0.1012865073234563);
01466     pnt[2] = Point(0.1012865073234563 , 0.7974269853530872);
01467     pnt[3] = Point(0.1012865073234563 , 0.1012865073234563);
01468     pnt[4] = Point(0.0597158717897698 , 0.4701420641051151);
01469     pnt[5] = Point(0.4701420641051151 , 0.0597158717897698);
01470     pnt[6] = Point(0.4701420641051151 , 0.4701420641051151);
01471     pnt[7] = Point(0.5357953464498992 , 0.2321023267750504);
01472     pnt[8] = Point(0.2321023267750504 , 0.5357953464498992);
01473     pnt[9] = Point(0.2321023267750504 , 0.2321023267750504);
01474     pnt[10] = Point(0.9410382782311209 , 0.0294808608844396);
01475     pnt[11] = Point(0.0294808608844396 , 0.9410382782311209);
01476     pnt[12] = Point(0.0294808608844396 , 0.0294808608844396);
01477     pnt[13] = Point(0.7384168123405100 , 0.2321023267750504);
01478     pnt[14] = Point(0.7384168123405100 , 0.0294808608844396);
01479     pnt[15] = Point(0.2321023267750504 , 0.7384168123405100);
01480     pnt[16] = Point(0.2321023267750504 , 0.0294808608844396);
01481     pnt[17] = Point(0.0294808608844396 , 0.7384168123405100);
01482     pnt[18] = Point(0.0294808608844396 , 0.2321023267750504);
01483     weight[0] = 0.0378610912003147;
01484     weight[1] = 0.0376204254131829;
01485     weight[2] = 0.0376204254131829;
01486     weight[3] = 0.0376204254131829;
01487     weight[4] = 0.0783573522441174;
01488     weight[5] = 0.0783573522441174;
01489     weight[6] = 0.0783573522441174;
01490     weight[7] = 0.1162714796569659;
01491     weight[8] = 0.1162714796569659;
01492     weight[9] = 0.1162714796569659;
01493     weight[10] = 0.0134442673751655;
01494     weight[11] = 0.0134442673751655;
01495     weight[12] = 0.0134442673751655;
01496     weight[13] = 0.0375097224552317;
01497     weight[14] = 0.0375097224552317;
01498     weight[15] = 0.0375097224552317;
01499     weight[16] = 0.0375097224552317;
01500     weight[17] = 0.0375097224552317;
01501     weight[18] = 0.0375097224552317;
01502     } else if (deg == 9) {
01503     pnt.resize(19); weight.resize(19);
01504     pnt[0] = Point(0.33333333333333331     ,  0.33333333333333331);
01505     pnt[1] = Point(2.06349616025259287E-002,  0.48968251919873701);
01506     pnt[2] = Point(0.48968251919873701     ,  2.06349616025259287E-002);
01507     pnt[3] = Point(0.48968251919873701      , 0.48968251919873701);
01508     pnt[4] = Point(0.12582081701412900     ,  0.43708959149293553);
01509     pnt[5] = Point(0.43708959149293553     ,  0.12582081701412900);
01510     pnt[6] = Point(0.43708959149293553     ,  0.43708959149293553);
01511     pnt[7] = Point(0.62359292876193562     ,  0.18820353561903219);
01512     pnt[8] = Point(0.18820353561903219     ,  0.62359292876193562);
01513     pnt[9] = Point(0.18820353561903219     ,  0.18820353561903219);
01514     pnt[10] = Point(0.91054097321109406     ,  4.47295133944529688E-002);
01515     pnt[11] = Point(4.47295133944529688E-002,  0.91054097321109406);
01516     pnt[12] = Point(4.47295133944529688E-002,  4.47295133944529688E-002);
01517     pnt[13] = Point(0.74119859878449801     ,  3.68384120547362581E-002);
01518     pnt[14] = Point(0.74119859878449801     ,  0.22196298916076573);
01519     pnt[15] = Point(3.68384120547362581E-002,  0.74119859878449801);
01520     pnt[16] = Point(3.68384120547362581E-002,  0.22196298916076573);
01521     pnt[17] = Point(0.22196298916076573     ,  0.74119859878449801);
01522     pnt[18] = Point(0.22196298916076573     ,  3.68384120547362581E-002 );
01523     weight[0] = 9.71357962827961025E-002;
01524     weight[1] = 3.13347002271398278E-002;
01525     weight[2] = 3.13347002271398278E-002;
01526     weight[3] = 3.13347002271398278E-002;
01527     weight[4] = 7.78275410047754301E-002;
01528     weight[5] = 7.78275410047754301E-002;
01529     weight[6] = 7.78275410047754301E-002;
01530     weight[7] = 7.96477389272090969E-002;
01531     weight[8] = 7.96477389272090969E-002;
01532     weight[9] = 7.96477389272090969E-002;
01533     weight[10] = 2.55776756586981006E-002;
01534     weight[11] = 2.55776756586981006E-002;
01535     weight[12] = 2.55776756586981006E-002;
01536     weight[13] = 4.32835393772893970E-002;
01537     weight[14] = 4.32835393772893970E-002;
01538     weight[15] = 4.32835393772893970E-002;
01539     weight[16] = 4.32835393772893970E-002;
01540     weight[17] = 4.32835393772893970E-002;
01541     weight[18] = 4.32835393772893970E-002;
01542     } else if (deg<=11) {
01543     pnt.resize(28); weight.resize(28);
01544     pnt[0] = Point(0.33333333333333333 , 0.333333333333333333);
01545     pnt[1] = Point(0.9480217181434233  , 0.02598914092828833);
01546     pnt[2] = Point(0.02598914092828833 , 0.9480217181434233);
01547     pnt[3] = Point(0.02598914092828833 , 0.02598914092828833);
01548     pnt[4] = Point(0.8114249947041546  , 0.09428750264792270);
01549     pnt[5] = Point(0.09428750264792270 , 0.8114249947041546);
01550     pnt[6] = Point(0.09428750264792270 , 0.09428750264792270);
01551     pnt[7] = Point(0.01072644996557060 , 0.4946367750172147);
01552     pnt[8] = Point(0.4946367750172147  , 0.01072644996557060);
01553     pnt[9] = Point(0.4946367750172147  , 0.4946367750172147);
01554     pnt[10] = Point(0.5853132347709715  , 0.2073433826145142);
01555     pnt[11] = Point(0.2073433826145142  , 0.5853132347709715);
01556     pnt[12] = Point(0.2073433826145142  , 0.2073433826145142);
01557     pnt[13] = Point(0.1221843885990187  , 0.4389078057004907);
01558     pnt[14] = Point(0.4389078057004907  , 0.1221843885990187);
01559     pnt[15] = Point(0.4389078057004907  , 0.4389078057004907);
01560     pnt[16] = Point(0.6779376548825902  , 0.04484167758913055);
01561     pnt[17] = Point(0.6779376548825902  , 0.27722066752827925);
01562     pnt[18] = Point(0.04484167758913055 , 0.6779376548825902);
01563     pnt[19] = Point(0.04484167758913055 , 0.27722066752827925);
01564     pnt[20] = Point(0.27722066752827925 , 0.6779376548825902);
01565     pnt[21] = Point(0.27722066752827925 , 0.04484167758913055);
01566     pnt[22] = Point(0.8588702812826364  , 0.00000000000000000);
01567     pnt[23] = Point(0.8588702812826364  , 0.1411297187173636);
01568     pnt[24] = Point(0.0000000000000000  , 0.8588702812826364);
01569     pnt[25] = Point(0.0000000000000000  , 0.1411297187173636);
01570     pnt[26] = Point(0.1411297187173636  , 0.8588702812826364);
01571     pnt[27] = Point(0.1411297187173636  , 0.0000000000000000);
01572     weight[0] = 0.08797730116222190;
01573     weight[1] = 0.008744311553736190;
01574     weight[2] = 0.008744311553736190;
01575     weight[3] = 0.008744311553736190;
01576     weight[4] = 0.03808157199393533;
01577     weight[5] = 0.03808157199393533;
01578     weight[6] = 0.03808157199393533;
01579     weight[7] = 0.01885544805613125;
01580     weight[8] = 0.01885544805613125;
01581     weight[9] = 0.01885544805613125;
01582     weight[10] = 0.07215969754474100;
01583     weight[11] = 0.07215969754474100;
01584     weight[12] = 0.07215969754474100;
01585     weight[13] = 0.06932913870553720;
01586     weight[14] = 0.06932913870553720;
01587     weight[15] = 0.06932913870553720;
01588     weight[16] = 0.04105631542928860;
01589     weight[17] = 0.04105631542928860;
01590     weight[18] = 0.04105631542928860;
01591     weight[19] = 0.04105631542928860;
01592     weight[20] = 0.04105631542928860;
01593     weight[21] = 0.04105631542928860;
01594     weight[22] = 0.007362383783300573;
01595     weight[23] = 0.007362383783300573;
01596     weight[24] = 0.007362383783300573;
01597     weight[25] = 0.007362383783300573;
01598     weight[26] = 0.007362383783300573;
01599     weight[27] = 0.007362383783300573;
01600     } else if (deg<=13) {
01601       pnt.resize(37); weight.resize(37);
01602       pnt[0] = Point(0.333333333333333333333333333333,  0.333333333333333333333333333333);
01603       pnt[1] = Point(0.950275662924105565450352089520,  0.024862168537947217274823955239);
01604       pnt[2] = Point(0.024862168537947217274823955239,  0.950275662924105565450352089520);
01605       pnt[3] = Point(0.024862168537947217274823955239,  0.024862168537947217274823955239);
01606       pnt[4] = Point(0.171614914923835347556304795551,  0.414192542538082326221847602214);
01607       pnt[5] = Point(0.414192542538082326221847602214,  0.171614914923835347556304795551);
01608       pnt[6] = Point(0.414192542538082326221847602214,  0.414192542538082326221847602214);
01609       pnt[7] = Point(0.539412243677190440263092985511,  0.230293878161404779868453507244);
01610       pnt[8] = Point(0.230293878161404779868453507244,  0.539412243677190440263092985511);
01611       pnt[9] = Point(0.230293878161404779868453507244,  0.230293878161404779868453507244);
01612       pnt[10] = Point(0.772160036676532561750285570113,  0.113919981661733719124857214943);
01613       pnt[11] = Point(0.113919981661733719124857214943,  0.772160036676532561750285570113);
01614       pnt[12] = Point(0.113919981661733719124857214943,  0.113919981661733719124857214943);
01615       pnt[13] = Point(0.009085399949835353883572964740,  0.495457300025082323058213517632);
01616       pnt[14] = Point(0.495457300025082323058213517632,  0.009085399949835353883572964740);
01617       pnt[15] = Point(0.495457300025082323058213517632,  0.495457300025082323058213517632);
01618       pnt[16] = Point(0.062277290305886993497083640527,  0.468861354847056503251458179727);
01619       pnt[17] = Point(0.468861354847056503251458179727,  0.062277290305886993497083640527);
01620       pnt[18] = Point(0.468861354847056503251458179727,  0.468861354847056503251458179727);
01621       pnt[19] = Point(0.022076289653624405142446876931,  0.851306504174348550389457672223);
01622       pnt[20] = Point(0.022076289653624405142446876931,  0.126617206172027096933163647918);
01623       pnt[21] = Point(0.851306504174348550389457672223,  0.022076289653624405142446876931);
01624       pnt[22] = Point(0.851306504174348550389457672223,  0.126617206172027096933163647918);
01625       pnt[23] = Point(0.126617206172027096933163647918,  0.022076289653624405142446876931);
01626       pnt[24] = Point(0.126617206172027096933163647918,  0.851306504174348550389457672223);
01627       pnt[25] = Point(0.018620522802520968955913511549,  0.689441970728591295496647976487);
01628       pnt[26] = Point(0.018620522802520968955913511549,  0.291937506468887771754472382212);
01629       pnt[27] = Point(0.689441970728591295496647976487,  0.018620522802520968955913511549);
01630       pnt[28] = Point(0.689441970728591295496647976487,  0.291937506468887771754472382212);
01631       pnt[29] = Point(0.291937506468887771754472382212,  0.018620522802520968955913511549);
01632       pnt[30] = Point(0.291937506468887771754472382212,  0.689441970728591295496647976487);
01633       pnt[31] = Point(0.096506481292159228736516560903,  0.635867859433872768286976979827);
01634       pnt[32] = Point(0.096506481292159228736516560903,  0.267625659273967961282458816185);
01635       pnt[33] = Point(0.635867859433872768286976979827,  0.096506481292159228736516560903);
01636       pnt[34] = Point(0.635867859433872768286976979827,  0.267625659273967961282458816185);
01637       pnt[35] = Point(0.267625659273967961282458816185,  0.096506481292159228736516560903);
01638       pnt[36] = Point(0.267625659273967961282458816185,  0.635867859433872768286976979827);
01639 
01640         weight[0] = 0.051739766065744133555179145422;
01641         weight[1] = 0.008007799555564801597804123460;
01642         weight[2] = 0.008007799555564801597804123460;
01643         weight[3] = 0.008007799555564801597804123460;
01644         weight[4] = 0.046868898981821644823226732071;
01645         weight[5] = 0.046868898981821644823226732071;
01646         weight[6] = 0.046868898981821644823226732071;
01647         weight[7] = 0.046590940183976487960361770070;
01648         weight[8] = 0.046590940183976487960361770070;
01649         weight[9] = 0.046590940183976487960361770070;
01650         weight[10] = 0.031016943313796381407646220131;
01651         weight[11] = 0.031016943313796381407646220131;
01652         weight[12] = 0.031016943313796381407646220131;
01653         weight[13] = 0.010791612736631273623178240136;
01654         weight[14] = 0.010791612736631273623178240136;
01655         weight[15] = 0.010791612736631273623178240136;
01656         weight[16] = 0.032195534242431618819414482205;
01657         weight[17] = 0.032195534242431618819414482205;
01658         weight[18] = 0.032195534242431618819414482205;
01659         weight[19] = 0.015445834210701583817692900053;
01660         weight[20] = 0.015445834210701583817692900053;
01661         weight[21] = 0.015445834210701583817692900053;
01662         weight[22] = 0.015445834210701583817692900053;
01663         weight[23] = 0.015445834210701583817692900053;
01664         weight[24] = 0.015445834210701583817692900053;
01665         weight[25] = 0.017822989923178661888748319485;
01666         weight[26] = 0.017822989923178661888748319485;
01667         weight[27] = 0.017822989923178661888748319485;
01668         weight[28] = 0.017822989923178661888748319485;
01669         weight[29] = 0.017822989923178661888748319485;
01670         weight[30] = 0.017822989923178661888748319485;
01671         weight[31] = 0.037038683681384627918546472190;
01672         weight[32] = 0.037038683681384627918546472190;
01673         weight[33] = 0.037038683681384627918546472190;
01674         weight[34] = 0.037038683681384627918546472190;
01675         weight[35] = 0.037038683681384627918546472190;
01676         weight[36] = 0.037038683681384627918546472190;
01677     } else if (deg<=14) {
01678       pnt.resize(46); weight.resize(46);
01679         pnt[0] = Point(0.33333333333333331000 , 0.33333333333333331000);
01680         pnt[1] = Point(0.00997976080645843200 , 0.00997976080645843200);
01681         pnt[2] = Point(0.00997976080645843200 , 0.98004047838708308000);
01682         pnt[3] = Point(0.98004047838708308000 , 0.00997976080645843200);
01683         pnt[4] = Point(0.47997789352118841000 , 0.47997789352118841000);
01684         pnt[5] = Point(0.47997789352118841000 , 0.04004421295762317100);
01685         pnt[6] = Point(0.04004421295762317100 , 0.47997789352118841000);
01686         pnt[7] = Point(0.15381195917696691000 , 0.15381195917696691000);
01687         pnt[8] = Point(0.15381195917696691000 , 0.69237608164606623000);
01688         pnt[9] = Point(0.69237608164606623000 , 0.15381195917696691000);
01689         pnt[10] = Point(0.07402347711698781300 , 0.07402347711698781300);
01690         pnt[11] = Point(0.07402347711698781300 , 0.85195304576602437000);
01691         pnt[12] = Point(0.85195304576602437000 , 0.07402347711698781300);
01692         pnt[13] = Point(0.13035468250332999000 , 0.13035468250332999000);
01693         pnt[14] = Point(0.13035468250332999000 , 0.73929063499334002000);
01694         pnt[15] = Point(0.73929063499334002000 , 0.13035468250332999000);
01695         pnt[16] = Point(0.23061722602665313000 , 0.23061722602665313000);
01696         pnt[17] = Point(0.23061722602665313000 , 0.53876554794669373000);
01697         pnt[18] = Point(0.53876554794669373000 , 0.23061722602665313000);
01698         pnt[19] = Point(0.42233208341914780000 , 0.42233208341914780000);
01699         pnt[20] = Point(0.42233208341914780000 , 0.15533583316170441000);
01700         pnt[21] = Point(0.15533583316170441000 , 0.42233208341914780000);
01701         pnt[22] = Point(0.78623738593466097000 , 0.19061636003190091000);
01702         pnt[23] = Point(0.78623738593466097000 , 0.02314625403343817400);
01703         pnt[24] = Point(0.19061636003190091000 , 0.78623738593466097000);
01704         pnt[25] = Point(0.19061636003190091000 , 0.02314625403343817400);
01705         pnt[26] = Point(0.02314625403343817400 , 0.78623738593466097000);
01706         pnt[27] = Point(0.02314625403343817400 , 0.19061636003190091000);
01707         pnt[28] = Point(0.63055214366060741000 , 0.36232313774354713000);
01708         pnt[29] = Point(0.63055214366060741000 , 0.00712471859584540290);
01709         pnt[30] = Point(0.36232313774354713000 , 0.63055214366060741000);
01710         pnt[31] = Point(0.36232313774354713000 , 0.00712471859584540290);
01711         pnt[32] = Point(0.00712471859584540290 , 0.63055214366060741000);
01712         pnt[33] = Point(0.00712471859584540290 , 0.36232313774354713000);
01713         pnt[34] = Point(0.62657732985630632000 , 0.29077120588366739000);
01714         pnt[35] = Point(0.62657732985630632000 , 0.08265146426002623100);
01715         pnt[36] = Point(0.29077120588366739000 , 0.62657732985630632000);
01716         pnt[37] = Point(0.29077120588366739000 , 0.08265146426002623100);
01717         pnt[38] = Point(0.08265146426002623100 , 0.62657732985630632000);
01718         pnt[39] = Point(0.08265146426002623100 , 0.29077120588366739000);
01719         pnt[40] = Point(0.91420998492962546000 , 0.07116571087775076800);
01720         pnt[41] = Point(0.91420998492962546000 , 0.01462430419262372700);
01721         pnt[42] = Point(0.07116571087775076800 , 0.91420998492962546000);
01722         pnt[43] = Point(0.07116571087775076800 , 0.01462430419262372700);
01723         pnt[44] = Point(0.01462430419262372700 , 0.91420998492962546000);
01724         pnt[45] = Point(0.01462430419262372700 , 0.07116571087775076800);
01725 
01726         weight[0] = 0.05859628522602859700;
01727         weight[1] = 0.00173515122972526760;
01728         weight[2] = 0.00173515122972526760;
01729         weight[3] = 0.00173515122972526760;
01730         weight[4] = 0.02616378255861452300;
01731         weight[5] = 0.02616378255861452300;
01732         weight[6] = 0.02616378255861452300;
01733         weight[7] = 0.00391972924240182890;
01734         weight[8] = 0.00391972924240182890;
01735         weight[9] = 0.00391972924240182890;
01736         weight[10] = 0.01224735975694086700;
01737         weight[11] = 0.01224735975694086700;
01738         weight[12] = 0.01224735975694086700;
01739         weight[13] = 0.02819962850325796000;
01740         weight[14] = 0.02819962850325796000;
01741         weight[15] = 0.02819962850325796000;
01742         weight[16] = 0.05088708718595948800;
01743         weight[17] = 0.05088708718595948800;
01744         weight[18] = 0.05088708718595948800;
01745         weight[19] = 0.05045343990160360000;
01746         weight[20] = 0.05045343990160360000;
01747         weight[21] = 0.05045343990160360000;
01748         weight[22] = 0.01706364421223345200;
01749         weight[23] = 0.01706364421223345200;
01750         weight[24] = 0.01706364421223345200;
01751         weight[25] = 0.01706364421223345200;
01752         weight[26] = 0.01706364421223345200;
01753         weight[27] = 0.01706364421223345200;
01754         weight[28] = 0.00968346642550660040;
01755         weight[29] = 0.00968346642550660040;
01756         weight[30] = 0.00968346642550660040;
01757         weight[31] = 0.00968346642550660040;
01758         weight[32] = 0.00968346642550660040;
01759         weight[33] = 0.00968346642550660040;
01760         weight[34] = 0.03638575592848500300;
01761         weight[35] = 0.03638575592848500300;
01762         weight[36] = 0.03638575592848500300;
01763         weight[37] = 0.03638575592848500300;
01764         weight[38] = 0.03638575592848500300;
01765         weight[39] = 0.03638575592848500300;
01766         weight[40] = 0.00696466337351841270;
01767         weight[41] = 0.00696466337351841270;
01768         weight[42] = 0.00696466337351841270;
01769         weight[43] = 0.00696466337351841270;
01770         weight[44] = 0.00696466337351841270;
01771         weight[45] = 0.00696466337351841270;
01772     } else if (deg<=20) {
01773       pnt.resize(88); weight.resize(88);
01774         pnt[0] = Point(0.33333333333333331000 , 0.33333333333333331000);
01775         pnt[1] = Point(0.21587430593299198000 , 0.21587430593299198000);
01776         pnt[2] = Point(0.21587430593299198000 , 0.56825138813401610000);
01777         pnt[3] = Point(0.56825138813401610000 , 0.21587430593299198000);
01778         pnt[4] = Point(0.07537676652974727200 , 0.07537676652974727200);
01779         pnt[5] = Point(0.07537676652974727200 , 0.84924646694050543000);
01780         pnt[6] = Point(0.84924646694050543000 , 0.07537676652974727200);
01781         pnt[7] = Point(0.01030082813722179300 , 0.01030082813722179300);
01782         pnt[8] = Point(0.01030082813722179300 , 0.97939834372555645000);
01783         pnt[9] = Point(0.97939834372555645000 , 0.01030082813722179300);
01784         pnt[10] = Point(0.49360221129870019000 , 0.49360221129870019000);
01785         pnt[11] = Point(0.49360221129870019000 , 0.01279557740259962300);
01786         pnt[12] = Point(0.01279557740259962300 , 0.49360221129870019000);
01787         pnt[13] = Point(0.46155093810692532000 , 0.46155093810692532000);
01788         pnt[14] = Point(0.46155093810692532000 , 0.07689812378614935300);
01789         pnt[15] = Point(0.07689812378614935300 , 0.46155093810692532000);
01790         pnt[16] = Point(0.32862140642423698000 , 0.42934057025821037000);
01791         pnt[17] = Point(0.32862140642423698000 , 0.24203802331755264000);
01792         pnt[18] = Point(0.42934057025821037000 , 0.32862140642423698000);
01793         pnt[19] = Point(0.42934057025821037000 , 0.24203802331755264000);
01794         pnt[20] = Point(0.24203802331755264000 , 0.32862140642423698000);
01795         pnt[21] = Point(0.24203802331755264000 , 0.42934057025821037000);
01796         pnt[22] = Point(0.26048036178656875000 , 0.10157753428096944000);
01797         pnt[23] = Point(0.26048036178656875000 , 0.63794210393246176000);
01798         pnt[24] = Point(0.10157753428096944000 , 0.26048036178656875000);
01799         pnt[25] = Point(0.10157753428096944000 , 0.63794210393246176000);
01800         pnt[26] = Point(0.63794210393246176000 , 0.26048036178656875000);
01801         pnt[27] = Point(0.63794210393246176000 , 0.10157753428096944000);
01802         pnt[28] = Point(0.13707423584645531000 , 0.71006597300113017000);
01803         pnt[29] = Point(0.13707423584645531000 , 0.15285979115241455000);
01804         pnt[30] = Point(0.71006597300113017000 , 0.13707423584645531000);
01805         pnt[31] = Point(0.71006597300113017000 , 0.15285979115241455000);
01806         pnt[32] = Point(0.15285979115241455000 , 0.13707423584645531000);
01807         pnt[33] = Point(0.15285979115241455000 , 0.71006597300113017000);
01808         pnt[34] = Point(0.14672694587229979000 , 0.49854547767841484000);
01809         pnt[35] = Point(0.14672694587229979000 , 0.35472757644928543000);
01810         pnt[36] = Point(0.49854547767841484000 , 0.14672694587229979000);
01811         pnt[37] = Point(0.49854547767841484000 , 0.35472757644928543000);
01812         pnt[38] = Point(0.35472757644928543000 , 0.14672694587229979000);
01813         pnt[39] = Point(0.35472757644928543000 , 0.49854547767841484000);
01814         pnt[40] = Point(0.02699897774255329000 , 0.04918672267258199900);
01815         pnt[41] = Point(0.02699897774255329000 , 0.92381429958486472000);
01816         pnt[42] = Point(0.04918672267258199900 , 0.02699897774255329000);
01817         pnt[43] = Point(0.04918672267258199900 , 0.92381429958486472000);
01818         pnt[44] = Point(0.92381429958486472000 , 0.02699897774255329000);
01819         pnt[45] = Point(0.92381429958486472000 , 0.04918672267258199900);
01820         pnt[46] = Point(0.06187178593361702900 , 0.77966014654056937000);
01821         pnt[47] = Point(0.06187178593361702900 , 0.15846806752581366000);
01822         pnt[48] = Point(0.77966014654056937000 , 0.06187178593361702900);
01823         pnt[49] = Point(0.77966014654056937000 , 0.15846806752581366000);
01824         pnt[50] = Point(0.15846806752581366000 , 0.06187178593361702900);
01825         pnt[51] = Point(0.15846806752581366000 , 0.77966014654056937000);
01826         pnt[52] = Point(0.04772436742762199700 , 0.37049153914954763000);
01827         pnt[53] = Point(0.04772436742762199700 , 0.58178409342283044000);
01828         pnt[54] = Point(0.37049153914954763000 , 0.04772436742762199700);
01829         pnt[55] = Point(0.37049153914954763000 , 0.58178409342283044000);
01830         pnt[56] = Point(0.58178409342283044000 , 0.04772436742762199700);
01831         pnt[57] = Point(0.58178409342283044000 , 0.37049153914954763000);
01832         pnt[58] = Point(0.12060051518636437000 , 0.86334694875475260000);
01833         pnt[59] = Point(0.12060051518636437000 , 0.01605253605888301600);
01834         pnt[60] = Point(0.86334694875475260000 , 0.12060051518636437000);
01835         pnt[61] = Point(0.86334694875475260000 , 0.01605253605888301600);
01836         pnt[62] = Point(0.01605253605888301600 , 0.12060051518636437000);
01837         pnt[63] = Point(0.01605253605888301600 , 0.86334694875475260000);
01838         pnt[64] = Point(0.00269714779670978760 , 0.05619493818774550000);
01839         pnt[65] = Point(0.00269714779670978760 , 0.94110791401554472000);
01840         pnt[66] = Point(0.05619493818774550000 , 0.00269714779670978760);
01841         pnt[67] = Point(0.05619493818774550000 , 0.94110791401554472000);
01842         pnt[68] = Point(0.94110791401554472000 , 0.00269714779670978760);
01843         pnt[69] = Point(0.94110791401554472000 , 0.05619493818774550000);
01844         pnt[70] = Point(0.00301563327794236250 , 0.20867500674842135000);
01845         pnt[71] = Point(0.00301563327794236250 , 0.78830935997363627000);
01846         pnt[72] = Point(0.20867500674842135000 , 0.00301563327794236250);
01847         pnt[73] = Point(0.20867500674842135000 , 0.78830935997363627000);
01848         pnt[74] = Point(0.78830935997363627000 , 0.00301563327794236250);
01849         pnt[75] = Point(0.78830935997363627000 , 0.20867500674842135000);
01850         pnt[76] = Point(0.02990537578845702000 , 0.72115124091203409000);
01851         pnt[77] = Point(0.02990537578845702000 , 0.24894338329950894000);
01852         pnt[78] = Point(0.72115124091203409000 , 0.02990537578845702000);
01853         pnt[79] = Point(0.72115124091203409000 , 0.24894338329950894000);
01854         pnt[80] = Point(0.24894338329950894000 , 0.02990537578845702000);
01855         pnt[81] = Point(0.24894338329950894000 , 0.72115124091203409000);
01856         pnt[82] = Point(0.00675665422246098880 , 0.64005544194054187000);
01857         pnt[83] = Point(0.00675665422246098880 , 0.35318790383699716000);
01858         pnt[84] = Point(0.64005544194054187000 , 0.00675665422246098880);
01859         pnt[85] = Point(0.64005544194054187000 , 0.35318790383699716000);
01860         pnt[86] = Point(0.35318790383699716000 , 0.00675665422246098880);
01861         pnt[87] = Point(0.35318790383699716000 , 0.64005544194054187000);
01862 
01863         weight[0] = 0.01253760799449665600;
01864         weight[1] = 0.02747186987642421400;
01865         weight[2] = 0.02747186987642421400;
01866         weight[3] = 0.02747186987642421400;
01867         weight[4] = 0.00976527227705142370;
01868         weight[5] = 0.00976527227705142370;
01869         weight[6] = 0.00976527227705142370;
01870         weight[7] = 0.00139841953539182350;
01871         weight[8] = 0.00139841953539182350;
01872         weight[9] = 0.00139841953539182350;
01873         weight[10] = 0.00929210262518518310;
01874         weight[11] = 0.00929210262518518310;
01875         weight[12] = 0.00929210262518518310;
01876         weight[13] = 0.01657787603236692700;
01877         weight[14] = 0.01657787603236692700;
01878         weight[15] = 0.01657787603236692700;
01879         weight[16] = 0.02066776234866507900;
01880         weight[17] = 0.02066776234866507900;
01881         weight[18] = 0.02066776234866507900;
01882         weight[19] = 0.02066776234866507900;
01883         weight[20] = 0.02066776234866507900;
01884         weight[21] = 0.02066776234866507900;
01885         weight[22] = 0.02082223552115450600;
01886         weight[23] = 0.02082223552115450600;
01887         weight[24] = 0.02082223552115450600;
01888         weight[25] = 0.02082223552115450600;
01889         weight[26] = 0.02082223552115450600;
01890         weight[27] = 0.02082223552115450600;
01891         weight[28] = 0.00956863841984906090;
01892         weight[29] = 0.00956863841984906090;
01893         weight[30] = 0.00956863841984906090;
01894         weight[31] = 0.00956863841984906090;
01895         weight[32] = 0.00956863841984906090;
01896         weight[33] = 0.00956863841984906090;
01897         weight[34] = 0.02445277096897246300;
01898         weight[35] = 0.02445277096897246300;
01899         weight[36] = 0.02445277096897246300;
01900         weight[37] = 0.02445277096897246300;
01901         weight[38] = 0.02445277096897246300;
01902         weight[39] = 0.02445277096897246300;
01903         weight[40] = 0.00315573063063053420;
01904         weight[41] = 0.00315573063063053420;
01905         weight[42] = 0.00315573063063053420;
01906         weight[43] = 0.00315573063063053420;
01907         weight[44] = 0.00315573063063053420;
01908         weight[45] = 0.00315573063063053420;
01909         weight[46] = 0.01213679636532129800;
01910         weight[47] = 0.01213679636532129800;
01911         weight[48] = 0.01213679636532129800;
01912         weight[49] = 0.01213679636532129800;
01913         weight[50] = 0.01213679636532129800;
01914         weight[51] = 0.01213679636532129800;
01915         weight[52] = 0.01496648014388644900;
01916         weight[53] = 0.01496648014388644900;
01917         weight[54] = 0.01496648014388644900;
01918         weight[55] = 0.01496648014388644900;
01919         weight[56] = 0.01496648014388644900;
01920         weight[57] = 0.01496648014388644900;
01921         weight[58] = 0.00632759332177773930;
01922         weight[59] = 0.00632759332177773930;
01923         weight[60] = 0.00632759332177773930;
01924         weight[61] = 0.00632759332177773930;
01925         weight[62] = 0.00632759332177773930;
01926         weight[63] = 0.00632759332177773930;
01927         weight[64] = 0.00134256031206369590;
01928         weight[65] = 0.00134256031206369590;
01929         weight[66] = 0.00134256031206369590;
01930         weight[67] = 0.00134256031206369590;
01931         weight[68] = 0.00134256031206369590;
01932         weight[69] = 0.00134256031206369590;
01933         weight[70] = 0.00277607691634755400;
01934         weight[71] = 0.00277607691634755400;
01935         weight[72] = 0.00277607691634755400;
01936         weight[73] = 0.00277607691634755400;
01937         weight[74] = 0.00277607691634755400;
01938         weight[75] = 0.00277607691634755400;
01939         weight[76] = 0.01073984447418494100;
01940         weight[77] = 0.01073984447418494100;
01941         weight[78] = 0.01073984447418494100;
01942         weight[79] = 0.01073984447418494100;
01943         weight[80] = 0.01073984447418494100;
01944         weight[81] = 0.01073984447418494100;
01945         weight[82] = 0.00536780573818745280;
01946         weight[83] = 0.00536780573818745280;
01947         weight[84] = 0.00536780573818745280;
01948         weight[85] = 0.00536780573818745280;
01949         weight[86] = 0.00536780573818745280;
01950         weight[87] = 0.00536780573818745280;
01951     } else if (deg<=25) {
01952       pnt.resize(126); weight.resize(126);
01953         pnt[0] = Point(0.02794648307316999900 , 0.48602675846340998000);
01954         pnt[1] = Point(0.48602675846340998000 , 0.48602675846340998000);
01955         pnt[2] = Point(0.48602675846340998000 , 0.02794648307316999900);
01956         pnt[3] = Point(0.13117860132765000000 , 0.43441069933616999000);
01957         pnt[4] = Point(0.43441069933616999000 , 0.43441069933616999000);
01958         pnt[5] = Point(0.43441069933616999000 , 0.13117860132765000000);
01959         pnt[6] = Point(0.22022172951207000000 , 0.38988913524396002000);
01960         pnt[7] = Point(0.38988913524396002000 , 0.38988913524396002000);
01961         pnt[8] = Point(0.38988913524396002000 , 0.22022172951207000000);
01962         pnt[9] = Point(0.40311353196039001000 , 0.29844323401980000000);
01963         pnt[10] = Point(0.29844323401980000000 , 0.29844323401980000000);
01964         pnt[11] = Point(0.29844323401980000000 , 0.40311353196039001000);
01965         pnt[12] = Point(0.53191165532525997000 , 0.23404417233736999000);
01966         pnt[13] = Point(0.23404417233736999000 , 0.23404417233736999000);
01967         pnt[14] = Point(0.23404417233736999000 , 0.53191165532525997000);
01968         pnt[15] = Point(0.69706333078196003000 , 0.15146833460902001000);
01969         pnt[16] = Point(0.15146833460902001000 , 0.15146833460902001000);
01970         pnt[17] = Point(0.15146833460902001000 , 0.69706333078196003000);
01971         pnt[18] = Point(0.77453221290801000000 , 0.11273389354599000000);
01972         pnt[19] = Point(0.11273389354599000000 , 0.11273389354599000000);
01973         pnt[20] = Point(0.11273389354599000000 , 0.77453221290801000000);
01974         pnt[21] = Point(0.84456861581694997000 , 0.07771569209152999500);
01975         pnt[22] = Point(0.07771569209152999500 , 0.07771569209152999500);
01976         pnt[23] = Point(0.07771569209152999500 , 0.84456861581694997000);
01977         pnt[24] = Point(0.93021381277141002000 , 0.03489309361430000000);
01978         pnt[25] = Point(0.03489309361430000000 , 0.03489309361430000000);
01979         pnt[26] = Point(0.03489309361430000000 , 0.93021381277141002000);
01980         pnt[27] = Point(0.98548363075812995000 , 0.00725818462093000010);
01981         pnt[28] = Point(0.00725818462093000010 , 0.00725818462093000010);
01982         pnt[29] = Point(0.00725818462093000010 , 0.98548363075812995000);
01983         pnt[30] = Point(0.00129235270443999990 , 0.22721445215336000000);
01984         pnt[31] = Point(0.22721445215336000000 , 0.77149319514218995000);
01985         pnt[32] = Point(0.77149319514218995000 , 0.00129235270443999990);
01986         pnt[33] = Point(0.22721445215336000000 , 0.00129235270443999990);
01987         pnt[34] = Point(0.77149319514218995000 , 0.22721445215336000000);
01988         pnt[35] = Point(0.00129235270443999990 , 0.77149319514218995000);
01989         pnt[36] = Point(0.00539970127212000000 , 0.43501055485356999000);
01990         pnt[37] = Point(0.43501055485356999000 , 0.55958974387431004000);
01991         pnt[38] = Point(0.55958974387431004000 , 0.00539970127212000000);
01992         pnt[39] = Point(0.43501055485356999000 , 0.00539970127212000000);
01993         pnt[40] = Point(0.55958974387431004000 , 0.43501055485356999000);
01994         pnt[41] = Point(0.00539970127212000000 , 0.55958974387431004000);
01995         pnt[42] = Point(0.00638400303398000030 , 0.32030959927219999000);
01996         pnt[43] = Point(0.32030959927219999000 , 0.67330639769381995000);
01997         pnt[44] = Point(0.67330639769381995000 , 0.00638400303398000030);
01998         pnt[45] = Point(0.32030959927219999000 , 0.00638400303398000030);
01999         pnt[46] = Point(0.67330639769381995000 , 0.32030959927219999000);
02000         pnt[47] = Point(0.00638400303398000030 , 0.67330639769381995000);
02001         pnt[48] = Point(0.00502821150199000020 , 0.09175032228000999700);
02002         pnt[49] = Point(0.09175032228000999700 , 0.90322146621800004000);
02003         pnt[50] = Point(0.90322146621800004000 , 0.00502821150199000020);
02004         pnt[51] = Point(0.09175032228000999700 , 0.00502821150199000020);
02005         pnt[52] = Point(0.90322146621800004000 , 0.09175032228000999700);
02006         pnt[53] = Point(0.00502821150199000020 , 0.90322146621800004000);
02007         pnt[54] = Point(0.00682675862178000040 , 0.03801083585871999800);
02008         pnt[55] = Point(0.03801083585871999800 , 0.95516240551949005000);
02009         pnt[56] = Point(0.95516240551949005000 , 0.00682675862178000040);
02010         pnt[57] = Point(0.03801083585871999800 , 0.00682675862178000040);
02011         pnt[58] = Point(0.95516240551949005000 , 0.03801083585871999800);
02012         pnt[59] = Point(0.00682675862178000040 , 0.95516240551949005000);
02013         pnt[60] = Point(0.01001619963993000000 , 0.15742521848530999000);
02014         pnt[61] = Point(0.15742521848530999000 , 0.83255858187475995000);
02015         pnt[62] = Point(0.83255858187475995000 , 0.01001619963993000000);
02016         pnt[63] = Point(0.15742521848530999000 , 0.01001619963993000000);
02017         pnt[64] = Point(0.83255858187475995000 , 0.15742521848530999000);
02018         pnt[65] = Point(0.01001619963993000000 , 0.83255858187475995000);
02019         pnt[66] = Point(0.02575781317339000100 , 0.23988965977853000000);
02020         pnt[67] = Point(0.23988965977853000000 , 0.73435252704807996000);
02021         pnt[68] = Point(0.73435252704807996000 , 0.02575781317339000100);
02022         pnt[69] = Point(0.23988965977853000000 , 0.02575781317339000100);
02023         pnt[70] = Point(0.73435252704807996000 , 0.23988965977853000000);
02024         pnt[71] = Point(0.02575781317339000100 , 0.73435252704807996000);
02025         pnt[72] = Point(0.03022789811992000100 , 0.36194311812606000000);
02026         pnt[73] = Point(0.36194311812606000000 , 0.60782898375401995000);
02027         pnt[74] = Point(0.60782898375401995000 , 0.03022789811992000100);
02028         pnt[75] = Point(0.36194311812606000000 , 0.03022789811992000100);
02029         pnt[76] = Point(0.60782898375401995000 , 0.36194311812606000000);
02030         pnt[77] = Point(0.03022789811992000100 , 0.60782898375401995000);
02031         pnt[78] = Point(0.03050499010715999900 , 0.08355196095483000100);
02032         pnt[79] = Point(0.08355196095483000100 , 0.88594304893801001000);
02033         pnt[80] = Point(0.88594304893801001000 , 0.03050499010715999900);
02034         pnt[81] = Point(0.08355196095483000100 , 0.03050499010715999900);
02035         pnt[82] = Point(0.88594304893801001000 , 0.08355196095483000100);
02036         pnt[83] = Point(0.03050499010715999900 , 0.88594304893801001000);
02037         pnt[84] = Point(0.04595654736257000200 , 0.14844322073242000000);
02038         pnt[85] = Point(0.14844322073242000000 , 0.80560023190500996000);
02039         pnt[86] = Point(0.80560023190500996000 , 0.04595654736257000200);
02040         pnt[87] = Point(0.14844322073242000000 , 0.04595654736257000200);
02041         pnt[88] = Point(0.80560023190500996000 , 0.14844322073242000000);
02042         pnt[89] = Point(0.04595654736257000200 , 0.80560023190500996000);
02043         pnt[90] = Point(0.06744280054027999800 , 0.28373970872753002000);
02044         pnt[91] = Point(0.28373970872753002000 , 0.64881749073218997000);
02045         pnt[92] = Point(0.64881749073218997000 , 0.06744280054027999800);
02046         pnt[93] = Point(0.28373970872753002000 , 0.06744280054027999800);
02047         pnt[94] = Point(0.64881749073218997000 , 0.28373970872753002000);
02048         pnt[95] = Point(0.06744280054027999800 , 0.64881749073218997000);
02049         pnt[96] = Point(0.07004509141591000500 , 0.40689937511878999000);
02050         pnt[97] = Point(0.40689937511878999000 , 0.52305553346529998000);
02051         pnt[98] = Point(0.52305553346529998000 , 0.07004509141591000500);
02052         pnt[99] = Point(0.40689937511878999000 , 0.07004509141591000500);
02053         pnt[100] = Point(0.52305553346529998000 , 0.40689937511878999000);
02054         pnt[101] = Point(0.07004509141591000500 , 0.52305553346529998000);
02055         pnt[102] = Point(0.08391152464012000000 , 0.19411398702488999000);
02056         pnt[103] = Point(0.19411398702488999000 , 0.72197448833499001000);
02057         pnt[104] = Point(0.72197448833499001000 , 0.08391152464012000000);
02058         pnt[105] = Point(0.19411398702488999000 , 0.08391152464012000000);
02059         pnt[106] = Point(0.72197448833499001000 , 0.19411398702488999000);
02060         pnt[107] = Point(0.08391152464012000000 , 0.72197448833499001000);
02061         pnt[108] = Point(0.12037553567714999000 , 0.32413434700069998000);
02062         pnt[109] = Point(0.32413434700069998000 , 0.55549011732214004000);
02063         pnt[110] = Point(0.55549011732214004000 , 0.12037553567714999000);
02064         pnt[111] = Point(0.32413434700069998000 , 0.12037553567714999000);
02065         pnt[112] = Point(0.55549011732214004000 , 0.32413434700069998000);
02066         pnt[113] = Point(0.12037553567714999000 , 0.55549011732214004000);
02067         pnt[114] = Point(0.14806689915737001000 , 0.22927748355597999000);
02068         pnt[115] = Point(0.22927748355597999000 , 0.62265561728664998000);
02069         pnt[116] = Point(0.62265561728664998000 , 0.14806689915737001000);
02070         pnt[117] = Point(0.22927748355597999000 , 0.14806689915737001000);
02071         pnt[118] = Point(0.62265561728664998000 , 0.22927748355597999000);
02072         pnt[119] = Point(0.14806689915737001000 , 0.62265561728664998000);
02073         pnt[120] = Point(0.19177186586733000000 , 0.32561812259598000000);
02074         pnt[121] = Point(0.32561812259598000000 , 0.48261001153668998000);
02075         pnt[122] = Point(0.48261001153668998000 , 0.19177186586733000000);
02076         pnt[123] = Point(0.32561812259598000000 , 0.19177186586733000000);
02077         pnt[124] = Point(0.48261001153668998000 , 0.32561812259598000000);
02078         pnt[125] = Point(0.19177186586733000000 , 0.48261001153668998000);
02079         weight[0] = 0.00800558188002041710;
02080         weight[1] = 0.00800558188002041710;
02081         weight[2] = 0.00800558188002041710;
02082         weight[3] = 0.01594707683239050100;
02083         weight[4] = 0.01594707683239050100;
02084         weight[5] = 0.01594707683239050100;
02085         weight[6] = 0.01310914123079553000;
02086         weight[7] = 0.01310914123079553000;
02087         weight[8] = 0.01310914123079553000;
02088         weight[9] = 0.01958300096563562000;
02089         weight[10] = 0.01958300096563562000;
02090         weight[11] = 0.01958300096563562000;
02091         weight[12] = 0.01647088544153727000;
02092         weight[13] = 0.01647088544153727000;
02093         weight[14] = 0.01647088544153727000;
02094         weight[15] = 0.00854727907409210020;
02095         weight[16] = 0.00854727907409210020;
02096         weight[17] = 0.00854727907409210020;
02097         weight[18] = 0.00816188585722649180;
02098         weight[19] = 0.00816188585722649180;
02099         weight[20] = 0.00816188585722649180;
02100         weight[21] = 0.00612114653998377910;
02101         weight[22] = 0.00612114653998377910;
02102         weight[23] = 0.00612114653998377910;
02103         weight[24] = 0.00290849826493666490;
02104         weight[25] = 0.00290849826493666490;
02105         weight[26] = 0.00290849826493666490;
02106         weight[27] = 0.00069227524566199629;
02107         weight[28] = 0.00069227524566199629;
02108         weight[29] = 0.00069227524566199629;
02109         weight[30] = 0.00124828919927739700;
02110         weight[31] = 0.00124828919927739700;
02111         weight[32] = 0.00124828919927739700;
02112         weight[33] = 0.00124828919927739700;
02113         weight[34] = 0.00124828919927739700;
02114         weight[35] = 0.00124828919927739700;
02115         weight[36] = 0.00340475290880302200;
02116         weight[37] = 0.00340475290880302200;
02117         weight[38] = 0.00340475290880302200;
02118         weight[39] = 0.00340475290880302200;
02119         weight[40] = 0.00340475290880302200;
02120         weight[41] = 0.00340475290880302200;
02121         weight[42] = 0.00335965432606405090;
02122         weight[43] = 0.00335965432606405090;
02123         weight[44] = 0.00335965432606405090;
02124         weight[45] = 0.00335965432606405090;
02125         weight[46] = 0.00335965432606405090;
02126         weight[47] = 0.00335965432606405090;
02127         weight[48] = 0.00171615653949675410;
02128         weight[49] = 0.00171615653949675410;
02129         weight[50] = 0.00171615653949675410;
02130         weight[51] = 0.00171615653949675410;
02131         weight[52] = 0.00171615653949675410;
02132         weight[53] = 0.00171615653949675410;
02133         weight[54] = 0.00148085631671560600;
02134         weight[55] = 0.00148085631671560600;
02135         weight[56] = 0.00148085631671560600;
02136         weight[57] = 0.00148085631671560600;
02137         weight[58] = 0.00148085631671560600;
02138         weight[59] = 0.00148085631671560600;
02139         weight[60] = 0.00351131261072868500;
02140         weight[61] = 0.00351131261072868500;
02141         weight[62] = 0.00351131261072868500;
02142         weight[63] = 0.00351131261072868500;
02143         weight[64] = 0.00351131261072868500;
02144         weight[65] = 0.00351131261072868500;
02145         weight[66] = 0.00739355014970648380;
02146         weight[67] = 0.00739355014970648380;
02147         weight[68] = 0.00739355014970648380;
02148         weight[69] = 0.00739355014970648380;
02149         weight[70] = 0.00739355014970648380;
02150         weight[71] = 0.00739355014970648380;
02151         weight[72] = 0.00798308747737655820;
02152         weight[73] = 0.00798308747737655820;
02153         weight[74] = 0.00798308747737655820;
02154         weight[75] = 0.00798308747737655820;
02155         weight[76] = 0.00798308747737655820;
02156         weight[77] = 0.00798308747737655820;
02157         weight[78] = 0.00435596261315804140;
02158         weight[79] = 0.00435596261315804140;
02159         weight[80] = 0.00435596261315804140;
02160         weight[81] = 0.00435596261315804140;
02161         weight[82] = 0.00435596261315804140;
02162         weight[83] = 0.00435596261315804140;
02163         weight[84] = 0.00736505670141783180;
02164         weight[85] = 0.00736505670141783180;
02165         weight[86] = 0.00736505670141783180;
02166         weight[87] = 0.00736505670141783180;
02167         weight[88] = 0.00736505670141783180;
02168         weight[89] = 0.00736505670141783180;
02169         weight[90] = 0.01096357284641955000;
02170         weight[91] = 0.01096357284641955000;
02171         weight[92] = 0.01096357284641955000;
02172         weight[93] = 0.01096357284641955000;
02173         weight[94] = 0.01096357284641955000;
02174         weight[95] = 0.01096357284641955000;
02175         weight[96] = 0.01174996174354112100;
02176         weight[97] = 0.01174996174354112100;
02177         weight[98] = 0.01174996174354112100;
02178         weight[99] = 0.01174996174354112100;
02179         weight[100] = 0.01174996174354112100;
02180         weight[101] = 0.01174996174354112100;
02181         weight[102] = 0.01001560071379857000;
02182         weight[103] = 0.01001560071379857000;
02183         weight[104] = 0.01001560071379857000;
02184         weight[105] = 0.01001560071379857000;
02185         weight[106] = 0.01001560071379857000;
02186         weight[107] = 0.01001560071379857000;
02187         weight[108] = 0.01330964078762868000;
02188         weight[109] = 0.01330964078762868000;
02189         weight[110] = 0.01330964078762868000;
02190         weight[111] = 0.01330964078762868000;
02191         weight[112] = 0.01330964078762868000;
02192         weight[113] = 0.01330964078762868000;
02193         weight[114] = 0.01415444650522614000;
02194         weight[115] = 0.01415444650522614000;
02195         weight[116] = 0.01415444650522614000;
02196         weight[117] = 0.01415444650522614000;
02197         weight[118] = 0.01415444650522614000;
02198         weight[119] = 0.01415444650522614000;
02199         weight[120] = 0.01488137956116801000;
02200         weight[121] = 0.01488137956116801000;
02201         weight[122] = 0.01488137956116801000;
02202         weight[123] = 0.01488137956116801000;
02203         weight[124] = 0.01488137956116801000;
02204         weight[125] = 0.01488137956116801000;
02205     } else if (deg<=30) {
02206       pnt.resize(175); weight.resize(175);
02207         pnt[0] = Point(0.33333333333332998000 , 0.33333333333332998000);
02208         pnt[1] = Point(0.00733011643276999980 , 0.49633494178361998000);
02209         pnt[2] = Point(0.49633494178361998000 , 0.49633494178361998000);
02210         pnt[3] = Point(0.49633494178361998000 , 0.00733011643276999980);
02211         pnt[4] = Point(0.08299567580295999500 , 0.45850216209852002000);
02212         pnt[5] = Point(0.45850216209852002000 , 0.45850216209852002000);
02213         pnt[6] = Point(0.45850216209852002000 , 0.08299567580295999500);
02214         pnt[7] = Point(0.15098095612540999000 , 0.42450952193729002000);
02215         pnt[8] = Point(0.42450952193729002000 , 0.42450952193729002000);
02216         pnt[9] = Point(0.42450952193729002000 , 0.15098095612540999000);
02217         pnt[10] = Point(0.23590585989217000000 , 0.38204707005392002000);
02218         pnt[11] = Point(0.38204707005392002000 , 0.38204707005392002000);
02219         pnt[12] = Point(0.38204707005392002000 , 0.23590585989217000000);
02220         pnt[13] = Point(0.43802430840785000000 , 0.28098784579607999000);
02221         pnt[14] = Point(0.28098784579607999000 , 0.28098784579607999000);
02222         pnt[15] = Point(0.28098784579607999000 , 0.43802430840785000000);
02223         pnt[16] = Point(0.54530204829192996000 , 0.22734897585402999000);
02224         pnt[17] = Point(0.22734897585402999000 , 0.22734897585402999000);
02225         pnt[18] = Point(0.22734897585402999000 , 0.54530204829192996000);
02226         pnt[19] = Point(0.65088177698254002000 , 0.17455911150872999000);
02227         pnt[20] = Point(0.17455911150872999000 , 0.17455911150872999000);
02228         pnt[21] = Point(0.17455911150872999000 , 0.65088177698254002000);
02229         pnt[22] = Point(0.75348314559713003000 , 0.12325842720143999000);
02230         pnt[23] = Point(0.12325842720143999000 , 0.12325842720143999000);
02231         pnt[24] = Point(0.12325842720143999000 , 0.75348314559713003000);
02232         pnt[25] = Point(0.83983154221560996000 , 0.08008422889220000200);
02233         pnt[26] = Point(0.08008422889220000200 , 0.08008422889220000200);
02234         pnt[27] = Point(0.08008422889220000200 , 0.83983154221560996000);
02235         pnt[28] = Point(0.90445106518420004000 , 0.04777446740790000000);
02236         pnt[29] = Point(0.04777446740790000000 , 0.04777446740790000000);
02237         pnt[30] = Point(0.04777446740790000000 , 0.90445106518420004000);
02238         pnt[31] = Point(0.95655897063971995000 , 0.02172051468014000000);
02239         pnt[32] = Point(0.02172051468014000000 , 0.02172051468014000000);
02240         pnt[33] = Point(0.02172051468014000000 , 0.95655897063971995000);
02241         pnt[34] = Point(0.99047064476913005000 , 0.00476467761544000030);
02242         pnt[35] = Point(0.00476467761544000030 , 0.00476467761544000030);
02243         pnt[36] = Point(0.00476467761544000030 , 0.99047064476913005000);
02244         pnt[37] = Point(0.00092537119334999999 , 0.41529527091330998000);
02245         pnt[38] = Point(0.41529527091330998000 , 0.58377935789334001000);
02246         pnt[39] = Point(0.58377935789334001000 , 0.00092537119334999999);
02247         pnt[40] = Point(0.41529527091330998000 , 0.00092537119334999999);
02248         pnt[41] = Point(0.58377935789334001000 , 0.41529527091330998000);
02249         pnt[42] = Point(0.00092537119334999999 , 0.58377935789334001000);
02250         pnt[43] = Point(0.00138592585556000010 , 0.06118990978535000100);
02251         pnt[44] = Point(0.06118990978535000100 , 0.93742416435909004000);
02252         pnt[45] = Point(0.93742416435909004000 , 0.00138592585556000010);
02253         pnt[46] = Point(0.06118990978535000100 , 0.00138592585556000010);
02254         pnt[47] = Point(0.93742416435909004000 , 0.06118990978535000100);
02255         pnt[48] = Point(0.00138592585556000010 , 0.93742416435909004000);
02256         pnt[49] = Point(0.00368241545590999990 , 0.16490869013691001000);
02257         pnt[50] = Point(0.16490869013691001000 , 0.83140889440718002000);
02258         pnt[51] = Point(0.83140889440718002000 , 0.00368241545590999990);
02259         pnt[52] = Point(0.16490869013691001000 , 0.00368241545590999990);
02260         pnt[53] = Point(0.83140889440718002000 , 0.16490869013691001000);
02261         pnt[54] = Point(0.00368241545590999990 , 0.83140889440718002000);
02262         pnt[55] = Point(0.00390322342416000000 , 0.02503506223199999900);
02263         pnt[56] = Point(0.02503506223199999900 , 0.97106171434384003000);
02264         pnt[57] = Point(0.97106171434384003000 , 0.00390322342416000000);
02265         pnt[58] = Point(0.02503506223199999900 , 0.00390322342416000000);
02266         pnt[59] = Point(0.97106171434384003000 , 0.02503506223199999900);
02267         pnt[60] = Point(0.00390322342416000000 , 0.97106171434384003000);
02268         pnt[61] = Point(0.00323324815500999980 , 0.30606446515109997000);
02269         pnt[62] = Point(0.30606446515109997000 , 0.69070228669389000000);
02270         pnt[63] = Point(0.69070228669389000000 , 0.00323324815500999980);
02271         pnt[64] = Point(0.30606446515109997000 , 0.00323324815500999980);
02272         pnt[65] = Point(0.69070228669389000000 , 0.30606446515109997000);
02273         pnt[66] = Point(0.00323324815500999980 , 0.69070228669389000000);
02274         pnt[67] = Point(0.00646743211223999980 , 0.10707328373022000000);
02275         pnt[68] = Point(0.10707328373022000000 , 0.88645928415754005000);
02276         pnt[69] = Point(0.88645928415754005000 , 0.00646743211223999980);
02277         pnt[70] = Point(0.10707328373022000000 , 0.00646743211223999980);
02278         pnt[71] = Point(0.88645928415754005000 , 0.10707328373022000000);
02279         pnt[72] = Point(0.00646743211223999980 , 0.88645928415754005000);
02280         pnt[73] = Point(0.00324747549132999980 , 0.22995754934557999000);
02281         pnt[74] = Point(0.22995754934557999000 , 0.76679497516308004000);
02282         pnt[75] = Point(0.76679497516308004000 , 0.00324747549132999980);
02283         pnt[76] = Point(0.22995754934557999000 , 0.00324747549132999980);
02284         pnt[77] = Point(0.76679497516308004000 , 0.22995754934557999000);
02285         pnt[78] = Point(0.00324747549132999980 , 0.76679497516308004000);
02286         pnt[79] = Point(0.00867509080675000000 , 0.33703663330577999000);
02287         pnt[80] = Point(0.33703663330577999000 , 0.65428827588745997000);
02288         pnt[81] = Point(0.65428827588745997000 , 0.00867509080675000000);
02289         pnt[82] = Point(0.33703663330577999000 , 0.00867509080675000000);
02290         pnt[83] = Point(0.65428827588745997000 , 0.33703663330577999000);
02291         pnt[84] = Point(0.00867509080675000000 , 0.65428827588745997000);
02292         pnt[85] = Point(0.01559702646731000000 , 0.05625657618206000200);
02293         pnt[86] = Point(0.05625657618206000200 , 0.92814639735062998000);
02294         pnt[87] = Point(0.92814639735062998000 , 0.01559702646731000000);
02295         pnt[88] = Point(0.05625657618206000200 , 0.01559702646731000000);
02296         pnt[89] = Point(0.92814639735062998000 , 0.05625657618206000200);
02297         pnt[90] = Point(0.01559702646731000000 , 0.92814639735062998000);
02298         pnt[91] = Point(0.01797672125369000100 , 0.40245137521239999000);
02299         pnt[92] = Point(0.40245137521239999000 , 0.57957190353390997000);
02300         pnt[93] = Point(0.57957190353390997000 , 0.01797672125369000100);
02301         pnt[94] = Point(0.40245137521239999000 , 0.01797672125369000100);
02302         pnt[95] = Point(0.57957190353390997000 , 0.40245137521239999000);
02303         pnt[96] = Point(0.01797672125369000100 , 0.57957190353390997000);
02304         pnt[97] = Point(0.01712424535389000000 , 0.24365470201083000000);
02305         pnt[98] = Point(0.24365470201083000000 , 0.73922105263528004000);
02306         pnt[99] = Point(0.73922105263528004000 , 0.01712424535389000000);
02307         pnt[100] = Point(0.24365470201083000000 , 0.01712424535389000000);
02308         pnt[101] = Point(0.73922105263528004000 , 0.24365470201083000000);
02309         pnt[102] = Point(0.01712424535389000000 , 0.73922105263528004000);
02310         pnt[103] = Point(0.02288340534658000000 , 0.16538958561452999000);
02311         pnt[104] = Point(0.16538958561452999000 , 0.81172700903887995000);
02312         pnt[105] = Point(0.81172700903887995000 , 0.02288340534658000000);
02313         pnt[106] = Point(0.16538958561452999000 , 0.02288340534658000000);
02314         pnt[107] = Point(0.81172700903887995000 , 0.16538958561452999000);
02315         pnt[108] = Point(0.02288340534658000000 , 0.81172700903887995000);
02316         pnt[109] = Point(0.03273759728777000200 , 0.09930187449584999800);
02317         pnt[110] = Point(0.09930187449584999800 , 0.86796052821639003000);
02318         pnt[111] = Point(0.86796052821639003000 , 0.03273759728777000200);
02319         pnt[112] = Point(0.09930187449584999800 , 0.03273759728777000200);
02320         pnt[113] = Point(0.86796052821639003000 , 0.09930187449584999800);
02321         pnt[114] = Point(0.03273759728777000200 , 0.86796052821639003000);
02322         pnt[115] = Point(0.03382101234234000100 , 0.30847833306904998000);
02323         pnt[116] = Point(0.30847833306904998000 , 0.65770065458860005000);
02324         pnt[117] = Point(0.65770065458860005000 , 0.03382101234234000100);
02325         pnt[118] = Point(0.30847833306904998000 , 0.03382101234234000100);
02326         pnt[119] = Point(0.65770065458860005000 , 0.30847833306904998000);
02327         pnt[120] = Point(0.03382101234234000100 , 0.65770065458860005000);
02328         pnt[121] = Point(0.03554761446001999900 , 0.46066831859210999000);
02329         pnt[122] = Point(0.46066831859210999000 , 0.50378406694787004000);
02330         pnt[123] = Point(0.50378406694787004000 , 0.03554761446001999900);
02331         pnt[124] = Point(0.46066831859210999000 , 0.03554761446001999900);
02332         pnt[125] = Point(0.50378406694787004000 , 0.46066831859210999000);
02333         pnt[126] = Point(0.03554761446001999900 , 0.50378406694787004000);
02334         pnt[127] = Point(0.05053979030687000300 , 0.21881529945393000000);
02335         pnt[128] = Point(0.21881529945393000000 , 0.73064491023919997000);
02336         pnt[129] = Point(0.73064491023919997000 , 0.05053979030687000300);
02337         pnt[130] = Point(0.21881529945393000000 , 0.05053979030687000300);
02338         pnt[131] = Point(0.73064491023919997000 , 0.21881529945393000000);
02339         pnt[132] = Point(0.05053979030687000300 , 0.73064491023919997000);
02340         pnt[133] = Point(0.05701471491573000000 , 0.37920955156026998000);
02341         pnt[134] = Point(0.37920955156026998000 , 0.56377573352399002000);
02342         pnt[135] = Point(0.56377573352399002000 , 0.05701471491573000000);
02343         pnt[136] = Point(0.37920955156026998000 , 0.05701471491573000000);
02344         pnt[137] = Point(0.56377573352399002000 , 0.37920955156026998000);
02345         pnt[138] = Point(0.05701471491573000000 , 0.56377573352399002000);
02346         pnt[139] = Point(0.06415280642119999800 , 0.14296081941819000000);
02347         pnt[140] = Point(0.14296081941819000000 , 0.79288637416061003000);
02348         pnt[141] = Point(0.79288637416061003000 , 0.06415280642119999800);
02349         pnt[142] = Point(0.14296081941819000000 , 0.06415280642119999800);
02350         pnt[143] = Point(0.79288637416061003000 , 0.14296081941819000000);
02351         pnt[144] = Point(0.06415280642119999800 , 0.79288637416061003000);
02352         pnt[145] = Point(0.08050114828762999800 , 0.28373128210592002000);
02353         pnt[146] = Point(0.28373128210592002000 , 0.63576756960644998000);
02354         pnt[147] = Point(0.63576756960644998000 , 0.08050114828762999800);
02355         pnt[148] = Point(0.28373128210592002000 , 0.08050114828762999800);
02356         pnt[149] = Point(0.63576756960644998000 , 0.28373128210592002000);
02357         pnt[150] = Point(0.08050114828762999800 , 0.63576756960644998000);
02358         pnt[151] = Point(0.10436706813453001000 , 0.19673744100443999000);
02359         pnt[152] = Point(0.19673744100443999000 , 0.69889549086102998000);
02360         pnt[153] = Point(0.69889549086102998000 , 0.10436706813453001000);
02361         pnt[154] = Point(0.19673744100443999000 , 0.10436706813453001000);
02362         pnt[155] = Point(0.69889549086102998000 , 0.19673744100443999000);
02363         pnt[156] = Point(0.10436706813453001000 , 0.69889549086102998000);
02364         pnt[157] = Point(0.11384489442875000000 , 0.35588914121165999000);
02365         pnt[158] = Point(0.35588914121165999000 , 0.53026596435958995000);
02366         pnt[159] = Point(0.53026596435958995000 , 0.11384489442875000000);
02367         pnt[160] = Point(0.35588914121165999000 , 0.11384489442875000000);
02368         pnt[161] = Point(0.53026596435958995000 , 0.35588914121165999000);
02369         pnt[162] = Point(0.11384489442875000000 , 0.53026596435958995000);
02370         pnt[163] = Point(0.14536348771551999000 , 0.25981868535190999000);
02371         pnt[164] = Point(0.25981868535190999000 , 0.59481782693256002000);
02372         pnt[165] = Point(0.59481782693256002000 , 0.14536348771551999000);
02373         pnt[166] = Point(0.25981868535190999000 , 0.14536348771551999000);
02374         pnt[167] = Point(0.59481782693256002000 , 0.25981868535190999000);
02375         pnt[168] = Point(0.14536348771551999000 , 0.59481782693256002000);
02376         pnt[169] = Point(0.18994565282198000000 , 0.32192318123129998000);
02377         pnt[170] = Point(0.32192318123129998000 , 0.48813116594672001000);
02378         pnt[171] = Point(0.48813116594672001000 , 0.18994565282198000000);
02379         pnt[172] = Point(0.32192318123129998000 , 0.18994565282198000000);
02380         pnt[173] = Point(0.48813116594672001000 , 0.32192318123129998000);
02381         pnt[174] = Point(0.18994565282198000000 , 0.48813116594672001000);
02382 
02383         weight[0] = 0.01557996020289920000;
02384         weight[1] = 0.00317723370053413400;
02385         weight[2] = 0.00317723370053413400;
02386         weight[3] = 0.00317723370053413400;
02387         weight[4] = 0.01048342663573077100;
02388         weight[5] = 0.01048342663573077100;
02389         weight[6] = 0.01048342663573077100;
02390         weight[7] = 0.01320945957774363000;
02391         weight[8] = 0.01320945957774363000;
02392         weight[9] = 0.01320945957774363000;
02393         weight[10] = 0.01497500696627149900;
02394         weight[11] = 0.01497500696627149900;
02395         weight[12] = 0.01497500696627149900;
02396         weight[13] = 0.01498790444338419000;
02397         weight[14] = 0.01498790444338419000;
02398         weight[15] = 0.01498790444338419000;
02399         weight[16] = 0.01333886474102166000;
02400         weight[17] = 0.01333886474102166000;
02401         weight[18] = 0.01333886474102166000;
02402         weight[19] = 0.01088917111390201100;
02403         weight[20] = 0.01088917111390201100;
02404         weight[21] = 0.01088917111390201100;
02405         weight[22] = 0.00818944066089346070;
02406         weight[23] = 0.00818944066089346070;
02407         weight[24] = 0.00818944066089346070;
02408         weight[25] = 0.00557538758860778510;
02409         weight[26] = 0.00557538758860778510;
02410         weight[27] = 0.00557538758860778510;
02411         weight[28] = 0.00319121647341197600;
02412         weight[29] = 0.00319121647341197600;
02413         weight[30] = 0.00319121647341197600;
02414         weight[31] = 0.00129671514432704500;
02415         weight[32] = 0.00129671514432704500;
02416         weight[33] = 0.00129671514432704500;
02417         weight[34] = 0.00029826282613491719;
02418         weight[35] = 0.00029826282613491719;
02419         weight[36] = 0.00029826282613491719;
02420         weight[37] = 0.00099890568507889641;
02421         weight[38] = 0.00099890568507889641;
02422         weight[39] = 0.00099890568507889641;
02423         weight[40] = 0.00099890568507889641;
02424         weight[41] = 0.00099890568507889641;
02425         weight[42] = 0.00099890568507889641;
02426         weight[43] = 0.00046285084917325331;
02427         weight[44] = 0.00046285084917325331;
02428         weight[45] = 0.00046285084917325331;
02429         weight[46] = 0.00046285084917325331;
02430         weight[47] = 0.00046285084917325331;
02431         weight[48] = 0.00046285084917325331;
02432         weight[49] = 0.00123445133638241290;
02433         weight[50] = 0.00123445133638241290;
02434         weight[51] = 0.00123445133638241290;
02435         weight[52] = 0.00123445133638241290;
02436         weight[53] = 0.00123445133638241290;
02437         weight[54] = 0.00123445133638241290;
02438         weight[55] = 0.00057071985224320615;
02439         weight[56] = 0.00057071985224320615;
02440         weight[57] = 0.00057071985224320615;
02441         weight[58] = 0.00057071985224320615;
02442         weight[59] = 0.00057071985224320615;
02443         weight[60] = 0.00057071985224320615;
02444         weight[61] = 0.00112694612587762410;
02445         weight[62] = 0.00112694612587762410;
02446         weight[63] = 0.00112694612587762410;
02447         weight[64] = 0.00112694612587762410;
02448         weight[65] = 0.00112694612587762410;
02449         weight[66] = 0.00112694612587762410;
02450         weight[67] = 0.00174786694940733710;
02451         weight[68] = 0.00174786694940733710;
02452         weight[69] = 0.00174786694940733710;
02453         weight[70] = 0.00174786694940733710;
02454         weight[71] = 0.00174786694940733710;
02455         weight[72] = 0.00174786694940733710;
02456         weight[73] = 0.00118281881503165690;
02457         weight[74] = 0.00118281881503165690;
02458         weight[75] = 0.00118281881503165690;
02459         weight[76] = 0.00118281881503165690;
02460         weight[77] = 0.00118281881503165690;
02461         weight[78] = 0.00118281881503165690;
02462         weight[79] = 0.00199083929467503380;
02463         weight[80] = 0.00199083929467503380;
02464         weight[81] = 0.00199083929467503380;
02465         weight[82] = 0.00199083929467503380;
02466         weight[83] = 0.00199083929467503380;
02467         weight[84] = 0.00199083929467503380;
02468         weight[85] = 0.00190041279503598000;
02469         weight[86] = 0.00190041279503598000;
02470         weight[87] = 0.00190041279503598000;
02471         weight[88] = 0.00190041279503598000;
02472         weight[89] = 0.00190041279503598000;
02473         weight[90] = 0.00190041279503598000;
02474         weight[91] = 0.00449836580881745120;
02475         weight[92] = 0.00449836580881745120;
02476         weight[93] = 0.00449836580881745120;
02477         weight[94] = 0.00449836580881745120;
02478         weight[95] = 0.00449836580881745120;
02479         weight[96] = 0.00449836580881745120;
02480         weight[97] = 0.00347871946027471890;
02481         weight[98] = 0.00347871946027471890;
02482         weight[99] = 0.00347871946027471890;
02483         weight[100] = 0.00347871946027471890;
02484         weight[101] = 0.00347871946027471890;
02485         weight[102] = 0.00347871946027471890;
02486         weight[103] = 0.00410239903672395340;
02487         weight[104] = 0.00410239903672395340;
02488         weight[105] = 0.00410239903672395340;
02489         weight[106] = 0.00410239903672395340;
02490         weight[107] = 0.00410239903672395340;
02491         weight[108] = 0.00410239903672395340;
02492         weight[109] = 0.00402176154974416210;
02493         weight[110] = 0.00402176154974416210;
02494         weight[111] = 0.00402176154974416210;
02495         weight[112] = 0.00402176154974416210;
02496         weight[113] = 0.00402176154974416210;
02497         weight[114] = 0.00402176154974416210;
02498         weight[115] = 0.00603316466079506590;
02499         weight[116] = 0.00603316466079506590;
02500         weight[117] = 0.00603316466079506590;
02501         weight[118] = 0.00603316466079506590;
02502         weight[119] = 0.00603316466079506590;
02503         weight[120] = 0.00603316466079506590;
02504         weight[121] = 0.00394629030212959810;
02505         weight[122] = 0.00394629030212959810;
02506         weight[123] = 0.00394629030212959810;
02507         weight[124] = 0.00394629030212959810;
02508         weight[125] = 0.00394629030212959810;
02509         weight[126] = 0.00394629030212959810;
02510         weight[127] = 0.00664404453768026840;
02511         weight[128] = 0.00664404453768026840;
02512         weight[129] = 0.00664404453768026840;
02513         weight[130] = 0.00664404453768026840;
02514         weight[131] = 0.00664404453768026840;
02515         weight[132] = 0.00664404453768026840;
02516         weight[133] = 0.00825430585607845810;
02517         weight[134] = 0.00825430585607845810;
02518         weight[135] = 0.00825430585607845810;
02519         weight[136] = 0.00825430585607845810;
02520         weight[137] = 0.00825430585607845810;
02521         weight[138] = 0.00825430585607845810;
02522         weight[139] = 0.00649605663340641070;
02523         weight[140] = 0.00649605663340641070;
02524         weight[141] = 0.00649605663340641070;
02525         weight[142] = 0.00649605663340641070;
02526         weight[143] = 0.00649605663340641070;
02527         weight[144] = 0.00649605663340641070;
02528         weight[145] = 0.00925277814414660230;
02529         weight[146] = 0.00925277814414660230;
02530         weight[147] = 0.00925277814414660230;
02531         weight[148] = 0.00925277814414660230;
02532         weight[149] = 0.00925277814414660230;
02533         weight[150] = 0.00925277814414660230;
02534         weight[151] = 0.00916492072629427990;
02535         weight[152] = 0.00916492072629427990;
02536         weight[153] = 0.00916492072629427990;
02537         weight[154] = 0.00916492072629427990;
02538         weight[155] = 0.00916492072629427990;
02539         weight[156] = 0.00916492072629427990;
02540         weight[157] = 0.01156952462809767100;
02541         weight[158] = 0.01156952462809767100;
02542         weight[159] = 0.01156952462809767100;
02543         weight[160] = 0.01156952462809767100;
02544         weight[161] = 0.01156952462809767100;
02545         weight[162] = 0.01156952462809767100;
02546         weight[163] = 0.01176111646760917000;
02547         weight[164] = 0.01176111646760917000;
02548         weight[165] = 0.01176111646760917000;
02549         weight[166] = 0.01176111646760917000;
02550         weight[167] = 0.01176111646760917000;
02551         weight[168] = 0.01176111646760917000;
02552         weight[169] = 0.01382470218216540000;
02553         weight[170] = 0.01382470218216540000;
02554         weight[171] = 0.01382470218216540000;
02555         weight[172] = 0.01382470218216540000;
02556         weight[173] = 0.01382470218216540000;
02557         weight[174] = 0.01382470218216540000;
02558     }
02559     else {
02560       // throw error , we do not have have the required order
02561       TEUCHOS_TEST_FOR_EXCEPTION( true , std::runtime_error  , " GaussLobattoQuadrature::getTriangleQuadPoints order of the quadrature to high !!! ");
02562     }
02563 }