NGSolve  4.9
fem/intrule.hpp
00001 #ifndef FILE_INTRULE
00002 #define FILE_INTRULE
00003 
00004 /*********************************************************************/
00005 /* File:   intrule.hpp                                               */
00006 /* Author: Joachim Schoeberl                                         */
00007 /* Date:   25. Mar. 2000                                             */
00008 /*********************************************************************/
00009 
00010 namespace ngfem
00011 {
00012 
00014   class IntegrationPoint
00015   {
00016   private:
00018     int nr; 
00020     double pi[3];
00022     double weight;
00024     int facetnr;
00025 
00026   public:
00028     bool precomputed_geometry;
00029   public:
00031     IntegrationPoint & operator=(const IntegrationPoint & aip)
00032     {
00033       nr = aip.Nr();
00034       pi[0] = aip(0);
00035       pi[1] = aip(1);
00036       pi[2] = aip(2);
00037       weight = aip.Weight();
00038       precomputed_geometry = aip.precomputed_geometry;
00039       facetnr = -1;
00040       return *this;
00041     }
00042 
00044     IntegrationPoint (const double api[3], double aw)
00045     {
00046       nr = -1;
00047       pi[0] = api[0];
00048       pi[1] = api[1];
00049       pi[2] = api[2];
00050       weight = aw;
00051       facetnr = -1;
00052       precomputed_geometry = 0;
00053     }
00054 
00056     IntegrationPoint (double p1 = 0, double p2 = 0, double p3 = 0, double aw = 0)
00057     {
00058       nr = -1;
00059       pi[0] = p1;
00060       pi[1] = p2;
00061       pi[2] = p3;
00062       weight = aw;
00063       facetnr = -1;
00064       precomputed_geometry = 0;
00065     }
00066 
00068     IntegrationPoint (const FlatVector<double> & ap, double aw)
00069     {
00070       nr = -1;
00071       pi[0] = (ap.Size() >= 1) ? ap(0) : 0;
00072       pi[1] = (ap.Size() >= 2) ? ap(1) : 0;
00073       pi[2] = (ap.Size() >= 3) ? ap(2) : 0;
00074       weight = aw;
00075       facetnr = -1;
00076       precomputed_geometry = 0;
00077     }
00078 
00079 
00080     template <int D>
00081     IntegrationPoint (const Vec<D> & ap, double aw = -1)
00082     {
00083       nr = -1;
00084       for (int j = 0; j < D; j++)
00085         pi[j] = ap(j);
00086       weight = aw;
00087       facetnr = -1;
00088       precomputed_geometry = 0;
00089     }
00090 
00092     IntegrationPoint (const IntegrationPoint & aip)
00093     { *this = aip; }
00094 
00096     void SetNr (int anr) { nr = anr; }
00098     const double * Point () const { return pi; }
00100     int Size() const { return 3; }
00102     const double & operator() (int i) const { return pi[i]; }
00104     double & operator() (int i) { return pi[i]; }
00106     double Weight () const { return weight; }
00108     void SetWeight (double w) { weight = w; }
00109 
00111     int Nr () const { return nr; }
00112 
00114     int IPNr () const { return -1; }
00115   public:
00116     int & FacetNr() { return facetnr; }
00117     int FacetNr() const { return facetnr; }
00118 
00120     friend NGS_DLL_HEADER ostream & operator<< (ostream & ost, const IntegrationPoint & ip);
00121   };
00122 
00123 
00124 
00125   class NGS_DLL_HEADER ElementTransformation;
00126 
00132   class BaseMappedIntegrationPoint
00133   {
00134   protected:
00136     const IntegrationPoint * ip;
00138     const ElementTransformation * eltrans;
00140     double measure; 
00141   public:
00143     BaseMappedIntegrationPoint (const IntegrationPoint & aip,
00144                                 const ElementTransformation & aeltrans)
00145       : ip(&aip), eltrans(&aeltrans)  { ; }
00147     const IntegrationPoint & IP () const { return *ip; }
00149     const ElementTransformation & GetTransformation () const { return *eltrans; }
00151     int GetIPNr() const { return ip->Nr(); }
00153     double GetMeasure() const { return measure; }
00155     double GetWeight() const { return measure * ip->Weight(); }
00156   };
00157 
00158 
00159   template <int R, typename SCAL = double>
00160   class DimMappedIntegrationPoint : public BaseMappedIntegrationPoint
00161   {
00162   protected:
00164     Vec<R,SCAL> point;
00165   public:
00167     DimMappedIntegrationPoint (const IntegrationPoint & aip,
00168                                  const ElementTransformation & aeltrans)
00169       : BaseMappedIntegrationPoint (aip, aeltrans)
00170     { ; }
00172     const Vec<R,SCAL> & GetPoint () const { return point; }
00173     Vec<R,SCAL> & Point () { return point; }
00175     SCAL operator() (int i) const { return point(i); }
00176   };
00177 
00178 
00180   template <int DIMS, int DIMR, typename SCAL = double> 
00181   class MappedIntegrationPoint : public DimMappedIntegrationPoint<DIMR,SCAL>
00182   {
00183   private:
00185     Mat<DIMR,DIMS,SCAL> dxdxi;
00187     Mat<DIMS,DIMR,SCAL> dxidx;
00189     SCAL det;
00191     Vec<DIMR,SCAL> normalvec;
00192     Vec<DIMR,SCAL> tangentialvec;  // for independent integrator
00193   
00194  
00195   public:
00196     typedef SCAL TSCAL;
00198     NGS_DLL_HEADER MappedIntegrationPoint (const IntegrationPoint & aip,
00199                                            const ElementTransformation & aeltrans);
00200 
00201     MappedIntegrationPoint (const IntegrationPoint & aip,
00202                             const ElementTransformation & aeltrans,
00203                             int /* dummy */)
00204       : DimMappedIntegrationPoint<DIMR,SCAL> (aip, aeltrans)
00205     { ; }
00206 
00208     MappedIntegrationPoint (const IntegrationPoint & aip,
00209                             const ElementTransformation & aeltrans,
00210                             const FlatVec<DIMR, SCAL> ax,
00211                             const Mat<DIMR, DIMS, SCAL> & adxdxi)
00212       : DimMappedIntegrationPoint<DIMR,SCAL> (aip, aeltrans)
00213     {
00214       this->point = ax;
00215       dxdxi = adxdxi;
00216       Compute();
00217     }
00218 
00219     void Compute ()
00220     {
00221       if (DIMS == DIMR)
00222         {
00223           det = Det (dxdxi);
00224           dxidx = Inv (dxdxi);
00225           normalvec = TSCAL(0.0);
00226           tangentialvec = TSCAL(0.0);
00227         }
00228       else
00229         {
00230           if (DIMR == 3)
00231             {
00232               normalvec = Cross (Vec<3,SCAL> (dxdxi.Col(0)),
00233                                  Vec<3,SCAL> (dxdxi.Col(1)));
00234               det = L2Norm (normalvec);
00235               normalvec /= det;
00236             }
00237           else if (DIMR == 2)
00238             {
00239               det = sqrt ( sqr (dxdxi(0,0)) + sqr (dxdxi(1,0)));
00240 
00241               normalvec(0) = -dxdxi(1,0) / det;
00242               normalvec(1) = dxdxi(0,0) / det;
00243             }
00244           else
00245             {
00246               det = 1.0;
00247               normalvec = 1.0;
00248             }
00249         
00250           Mat<DIMS,DIMS,SCAL> ata, iata;
00251         
00252           ata = Trans (dxdxi) * dxdxi;
00253           iata = Inv (ata);
00254           dxidx = iata * Trans (dxdxi);
00255           tangentialvec = TSCAL(0.0);
00256         }
00257       this->measure = fabs (det);
00258     }
00259   
00261     const Mat<DIMR,DIMS,SCAL> & GetJacobian() const { return dxdxi; }
00262     Mat<DIMR,DIMS,SCAL> & Jacobian() { return dxdxi; }
00264     SCAL GetJacobiDet() const { return det; }
00266     const Mat<DIMS,DIMR,SCAL> & GetJacobianInverse() const { return dxidx; }
00268     const Vec<DIMR,SCAL> GetNV () const { return normalvec; }
00270     void SetNV ( const Vec<DIMR,SCAL> & vec) { normalvec = vec; }
00272     void SetTV ( const Vec<DIMR,SCAL> & vec) { tangentialvec = vec; }
00274     const Vec<DIMR,SCAL> GetTV () const { return tangentialvec; }
00276     int IsBoundary () const { return DIMS != DIMR; }
00277 
00279     void CalcHesse (Mat<2> & ddx1, Mat<2> & ddx2) const;
00280     void CalcHesse (Mat<2> & ddx1, Mat<2> & ddx2, Mat<2> & ddx3) const; 
00281     void CalcHesse (Mat<3> & ddx1, Mat<3> & ddx2, Mat<3> & ddx3) const;
00282   };
00283 
00284 
00285   template <int DIMS, int DIMR, typename SCAL> 
00286   inline ostream & operator<< (ostream & ost, const MappedIntegrationPoint<DIMS,DIMR,SCAL> & sip)
00287   {
00288     ost << sip.GetPoint() << ", dxdxi = " << sip.GetJacobian();
00289     return ost;
00290   }
00291 
00292 
00293 
00294 
00295 
00301   class IntegrationRule : public Array<IntegrationPoint> 
00302   {
00303   public:
00305     IntegrationRule () { ; }
00306 
00311     NGS_DLL_HEADER IntegrationRule (ELEMENT_TYPE eltype, int order);
00312 
00313 
00314     NGS_DLL_HEADER IntegrationRule (const IntegrationRule & ir2)
00315     {
00316       throw Exception ("should not call copy-constructor of ir\n");
00317     }
00318 
00319     NGS_DLL_HEADER IntegrationRule (int asize, LocalHeap & lh)
00320       : Array<IntegrationPoint> (asize, lh)
00321     {
00322       ;
00323     }
00324 
00325     NGS_DLL_HEADER IntegrationRule (int asize, double (*pts)[3], double * weights);
00326 
00327     // make it polymorphic
00328     virtual ~IntegrationRule() { ; }
00329 
00331     void AddIntegrationPoint (const IntegrationPoint & ip)
00332     { 
00333       Append (ip);
00334     }
00335 
00337     int GetNIP() const { return Size(); }
00338   };
00339 
00340   inline ostream & operator<< (ostream & ost, const IntegrationRule & ir)
00341   {
00342     for (int i = 0; i < ir.GetNIP(); i++)
00343       ost << ir[i] << endl;
00344     return ost;
00345   }
00346 
00347 
00348 
00349 
00350   template <int D>
00351   class IntegrationRuleTP : public IntegrationRule
00352   {
00353     const IntegrationRule *irx, *iry, *irz;
00354     // int sort[8];
00355 
00356     // ArrayMem<Vec<D>, 100> xi;
00357     // ArrayMem<double, 100> weight;
00358     ArrayMem<Vec<D>, 100> x;
00359     ArrayMem<Mat<D,D>, 100> dxdxi;
00360     ArrayMem<Mat<D,D>, 100> dxdxi_duffy;
00361 
00362   public:
00363     NGS_DLL_HEADER IntegrationRuleTP (const ElementTransformation & eltrans,
00364                        int order, bool compute_mapping, LocalHeap & lh);
00365 
00366     // tensor product rule for a facet
00367     NGS_DLL_HEADER IntegrationRuleTP (ELEMENT_TYPE eltype, FlatArray<int> sort, 
00368                        NODE_TYPE nt, int nodenr, int order, LocalHeap & lh);
00369 
00370     const IntegrationRule & GetIRX() const { return *irx; }
00371     const IntegrationRule & GetIRY() const { return *iry; }
00372     const IntegrationRule & GetIRZ() const { return *irz; }
00373 
00374     int GetNIP() const { return Size(); }
00375     double GetWeight (int i) const { return (*this)[i].Weight(); }  // weight[i]; }
00376     const Vec<D> GetXi(int i) const 
00377     { 
00378       Vec<D> ret;
00379       for (int j = 0; j < D; j++)
00380         ret(j) = (*this)[i](j);
00381       return ret;  // return xi[i]; 
00382     }
00383     Vec<D> & GetPoint(int i) const { return x[i]; }
00384     Mat<D,D> & GetJacobian(int i) const { return dxdxi[i]; }
00385     Mat<D,D> & GetDuffyJacobian(int i) const { return dxdxi_duffy[i]; }
00386   };
00387 
00388 
00389 
00390 
00396   NGS_DLL_HEADER extern void ComputeGaussRule (int n, 
00397                                 Array<double> & xi, 
00398                                 Array<double> & wi);
00399 
00400 
00401 
00402 
00409   NGS_DLL_HEADER extern void ComputeGaussJacobiRule (int n, 
00410                                       Array<double> & xi, 
00411                                       Array<double> & wi,
00412                                       double alf,
00413                                       double bet);
00414 
00415 
00421   NGS_DLL_HEADER extern void ComputeHermiteRule (int n, 
00422                                   Array<double> & x,
00423                                   Array<double> & w);
00424   
00425 
00426   extern NGS_DLL_HEADER const IntegrationRule & SelectIntegrationRule (ELEMENT_TYPE eltype, int order);
00427   extern NGS_DLL_HEADER const IntegrationRule & SelectIntegrationRuleJacobi10 (int order);
00428   extern NGS_DLL_HEADER const IntegrationRule & SelectIntegrationRuleJacobi20 (int order);
00429 
00430 
00431   // transformation of (d-1) dimensional integration points on facets to 
00432   // d-dimensional point in volumes
00433   class Facet2ElementTrafo
00434   {
00435   protected:
00436     // mutable IntegrationPoint elpoint;  
00437     ELEMENT_TYPE eltype;
00438     // const POINT3D * points;
00439     FlatVector<Vec<3> > points;
00440     const EDGE * edges;
00441     const FACE * faces;
00442     EDGE hedges[4];
00443     FACE hfaces[6];
00444   public:
00445     Facet2ElementTrafo(ELEMENT_TYPE aeltype) 
00446       : eltype(aeltype),
00447         points(99,(double*)ElementTopology::GetVertices (aeltype))
00448     {
00449       // points = ElementTopology::GetVertices (eltype);
00450       edges = ElementTopology::GetEdges (eltype);
00451       faces = ElementTopology::GetFaces (eltype);
00452     }
00453   
00454     Facet2ElementTrafo(ELEMENT_TYPE aeltype, const FlatArray<int> & vnums) 
00455       : eltype(aeltype),
00456         points(99,(double*)ElementTopology::GetVertices (aeltype))
00457     {
00458       // points = ElementTopology::GetVertices (eltype);
00459       edges = ElementTopology::GetEdges (eltype);
00460       faces = ElementTopology::GetFaces (eltype);
00461 
00462       if (eltype == ET_TRIG)
00463         {
00464           for (int i = 0; i < 3; i++)
00465             {
00466               hedges[i][0] = edges[i][0];
00467               hedges[i][1] = edges[i][1];
00468               if (vnums[hedges[i][0]] > vnums[hedges[i][1]])
00469                 swap (hedges[i][0], hedges[i][1]);
00470             }
00471           edges = &hedges[0];
00472         }
00473 
00474       if (eltype == ET_QUAD)
00475         {
00476           for (int i = 0; i < 4; i++)
00477             {
00478               hedges[i][0] = edges[i][0];
00479               hedges[i][1] = edges[i][1];
00480               if (vnums[hedges[i][0]] > vnums[hedges[i][1]])
00481                 swap (hedges[i][0], hedges[i][1]);
00482             }
00483           edges = &hedges[0];
00484         }
00485 
00486       if (eltype == ET_TET)
00487         {
00488           for (int i = 0; i < 4; i++)
00489             {
00490               hfaces[i][0] = faces[i][0];
00491               hfaces[i][1] = faces[i][1];
00492               hfaces[i][2] = faces[i][2];
00493               if (vnums[hfaces[i][0]] > vnums[hfaces[i][1]]) swap (hfaces[i][0], hfaces[i][1]);
00494               if (vnums[hfaces[i][1]] > vnums[hfaces[i][2]]) swap (hfaces[i][1], hfaces[i][2]);
00495               if (vnums[hfaces[i][0]] > vnums[hfaces[i][1]]) swap (hfaces[i][0], hfaces[i][1]);
00496             }
00497           faces = &hfaces[0];
00498         }
00499 
00500 
00501       if (eltype == ET_PRISM)
00502         {
00503           for (int i = 0; i < 2; i++)
00504             {
00505               hfaces[i][0] = faces[i][0];
00506               hfaces[i][1] = faces[i][1];
00507               hfaces[i][2] = faces[i][2];
00508               if (vnums[hfaces[i][0]] > vnums[hfaces[i][1]]) swap (hfaces[i][0], hfaces[i][1]);
00509               if (vnums[hfaces[i][1]] > vnums[hfaces[i][2]]) swap (hfaces[i][1], hfaces[i][2]);
00510               if (vnums[hfaces[i][0]] > vnums[hfaces[i][1]]) swap (hfaces[i][0], hfaces[i][1]);
00511             }
00512           for (int i = 2; i < 5; i++)
00513             {
00514               int jmin = 0;
00515               for (int j = 1; j < 4; j++)
00516                 if (vnums[faces[i][j]] < vnums[faces[i][jmin]]) jmin = j;
00517               int j1 = (jmin+1)%4;
00518               int j2 = (jmin+2)%4;
00519               int j3 = (jmin+3)%4;
00520               if (vnums[faces[i][j3]] < vnums[faces[i][j1]]) swap (j1, j3);
00521 
00522               hfaces[i][0] = faces[i][jmin];
00523               hfaces[i][1] = faces[i][j1];
00524               hfaces[i][2] = faces[i][j2];
00525               hfaces[i][3] = faces[i][j3];
00526             }
00527           faces = &hfaces[0];
00528         }
00529 
00530     }
00531 
00532 
00533     void operator()(int fnr, const IntegrationPoint &ipfac, IntegrationPoint & ipvol) const 
00534     {
00535       ELEMENT_TYPE facettype = ElementTopology::GetFacetType(eltype, fnr);
00536 
00537       switch (facettype)
00538         {
00539         case ET_SEGM:
00540           {
00541             FlatVec<3> p1 = points (edges[fnr][0]);
00542             FlatVec<3> p2 = points (edges[fnr][1]);
00543 
00544             // Vec<3> volp = p2 + ipfac(0) * (p1-p2);
00545             // ipvol = volp;
00546             ipvol = Vec<3> (p2 + ipfac(0) * (p1-p2));
00547             break;
00548           }
00549         case ET_TRIG:
00550           {
00551             FlatVec<3> p0 = points(faces[fnr][0]);
00552             FlatVec<3> p1 = points(faces[fnr][1]);
00553             FlatVec<3> p2 = points(faces[fnr][2]);
00554 
00555             // for (int j = 0; j < 3; j++)
00556             //   ipvol(j) = p2[j] + ipfac(0)*(p0[j]-p2[j]) + ipfac(1)*(p1[j]-p2[j]);
00557             ipvol = Vec<3> (p2 + ipfac(0) * (p0-p2) + ipfac(1)*(p1-p2));
00558             break;
00559           }
00560         case ET_QUAD:
00561           {
00562             FlatVec<3> p0 = points(faces[fnr][0]);
00563             FlatVec<3> p1 = points(faces[fnr][1]);
00564             FlatVec<3> p2 = points(faces[fnr][3]);
00565 
00566             // for (int j = 0; j < 3; j++)
00567             //  ipvol(j) = p0[j] + ipfac(0)*(p1[j]-p0[j]) + ipfac(1)*(p2[j]-p0[j]);
00568             ipvol = Vec<3> (p0 + ipfac(0) * (p1-p0) + ipfac(1)*(p2-p0));
00569             break;
00570           }
00571         default:
00572           throw Exception ("undefined facet type in Facet2ElementTrafo()\n");
00573         } 
00574       ipvol.FacetNr() = fnr;
00575     }
00576 
00577     const IntegrationPoint operator()(int fnr, const IntegrationPoint &ip1d) const 
00578     {
00579       IntegrationPoint elpoint;
00580       operator()(fnr, ip1d, elpoint);
00581       return elpoint;
00582     }
00583 
00584     IntegrationRule & operator() (int fnr, const IntegrationRule & irfacet, LocalHeap & lh)
00585     {
00586       IntegrationRule & irvol = *new (lh) IntegrationRule (irfacet.GetNIP(), lh);
00587 
00588       switch (ElementTopology::GetFacetType(eltype, fnr))
00589         {
00590         case ET_SEGM:
00591           {
00592             FlatVec<3> p1 = points (edges[fnr][0]);
00593             FlatVec<3> p2 = points (edges[fnr][1]);
00594 
00595             for (int i = 0; i < irfacet.GetNIP(); i++)
00596               irvol[i] = Vec<3> (p2 + irfacet[i](0) * (p1-p2));
00597             break;
00598           }
00599         case ET_TRIG:
00600           {
00601             FlatVec<3> p0 = points(faces[fnr][0]);
00602             FlatVec<3> p1 = points(faces[fnr][1]);
00603             FlatVec<3> p2 = points(faces[fnr][2]);
00604 
00605             for (int i = 0; i < irfacet.GetNIP(); i++)
00606               irvol[i] = Vec<3> (p2 + irfacet[i](0) * (p0-p2) + irfacet[i](1)*(p1-p2));
00607             break;
00608           }
00609         case ET_QUAD:
00610           {
00611             FlatVec<3> p0 = points(faces[fnr][0]);
00612             FlatVec<3> p1 = points(faces[fnr][1]);
00613             FlatVec<3> p2 = points(faces[fnr][3]);
00614 
00615             for (int i = 0; i < irfacet.GetNIP(); i++)
00616               irvol[i] = Vec<3> (p0 + irfacet[i](0) * (p1-p0) + irfacet[i](1)*(p2-p0));
00617             break;
00618           }
00619         default:
00620           throw Exception ("undefined facet type in Facet2ElementTrafo()\n");
00621         } 
00622 
00623       /*
00624       for (int i = 0; i < irfacet.GetNIP(); i++)
00625         (*this) (fnr, irfacet[i], irvol[i]);
00626       */
00627       for (int i = 0; i < irfacet.GetNIP(); i++)
00628         irvol[i].FacetNr() = fnr;
00629 
00630       return irvol;
00631     }
00632   };
00633 
00634 
00635 
00636 
00637 
00638 
00639 
00640   // transformation of (d-1) dimensional integration points on facets to 
00641   // d-dimensional point in volumes
00642   class Facet2SurfaceElementTrafo
00643   {
00644   protected:
00645     // mutable IntegrationPoint elpoint;  
00646     ELEMENT_TYPE eltype;
00647     const POINT3D * points;
00648     const EDGE * edges;
00649     const FACE * faces;
00650     EDGE hedges[4];
00651     FACE hfaces[6];
00652   public:
00653     Facet2SurfaceElementTrafo(ELEMENT_TYPE aeltype) : eltype(aeltype) 
00654     {
00655       points = ElementTopology::GetVertices (eltype);
00656       edges = ElementTopology::GetEdges (eltype);
00657       faces = ElementTopology::GetFaces (eltype);
00658     }
00659     
00660     Facet2SurfaceElementTrafo(ELEMENT_TYPE aeltype, FlatArray<int> & vnums) : eltype(aeltype) 
00661     {
00662       points = ElementTopology::GetVertices (eltype);
00663       edges = ElementTopology::GetEdges (eltype);
00664       faces = ElementTopology::GetFaces (eltype);
00665 
00666       if (eltype == ET_SEGM)
00667         {
00668           hedges[0][0] = edges[0][0];
00669           hedges[0][1] = edges[0][1];
00670           if (vnums[hedges[0][0]] > vnums[hedges[0][1]])
00671             swap (hedges[0][0], hedges[0][1]);
00672           edges = &hedges[0];
00673         }
00674 
00675       if (eltype == ET_TRIG)
00676         {
00677           hfaces[0][0] = faces[0][0];
00678           hfaces[0][1] = faces[0][1];
00679           hfaces[0][2] = faces[0][2];
00680           if (vnums[hfaces[0][0]] > vnums[hfaces[0][1]]) swap (hfaces[0][0], hfaces[0][1]);
00681           if (vnums[hfaces[0][1]] > vnums[hfaces[0][2]]) swap (hfaces[0][1], hfaces[0][2]);
00682           if (vnums[hfaces[0][0]] > vnums[hfaces[0][1]]) swap (hfaces[0][0], hfaces[0][1]);
00683           
00684           faces = &hfaces[0];
00685         }
00686 
00687       if (eltype == ET_QUAD)
00688         {
00689           int jmin = 0;
00690           for (int j = 1; j < 4; j++)
00691             if (vnums[faces[0][j]] < vnums[faces[0][jmin]]) jmin = j;
00692           int j1 = (jmin+1)%4;
00693           int j2 = (jmin+2)%4;
00694           int j3 = (jmin+3)%4;
00695           if (vnums[faces[0][j3]] < vnums[faces[0][j1]]) swap (j1, j3);
00696           
00697           hfaces[0][0] = faces[0][jmin];
00698           hfaces[0][1] = faces[0][j1];
00699           hfaces[0][2] = faces[0][j2];
00700           hfaces[0][3] = faces[0][j3];
00701           faces = &hfaces[0];
00702         }
00703     }
00704 
00705 
00706     void operator()(const IntegrationPoint &ipfac, IntegrationPoint & ipvol) const 
00707     {
00708       int fnr = 0;
00709       switch (eltype)
00710         {
00711         case ET_SEGM:
00712           {
00713             const POINT3D & p1 = points[edges[fnr][0]];
00714             const POINT3D & p2 = points[edges[fnr][1]];
00715             for (int j = 0; j < 3; j++)
00716               ipvol(j) = p2[j] + (ipfac(0))*(p1[j]-p2[j]);
00717             break;
00718           }
00719         case ET_TRIG:
00720           {
00721             const POINT3D & p0 = points[faces[fnr][0]];
00722             const POINT3D & p1 = points[faces[fnr][1]];
00723             const POINT3D & p2 = points[faces[fnr][2]];
00724             ipvol(0) = p2[0] + ipfac(0)*(p0[0]-p2[0]) + ipfac(1)*(p1[0]-p2[0]);
00725             ipvol(1) = p2[1] + ipfac(0)*(p0[1]-p2[1]) + ipfac(1)*(p1[1]-p2[1]);
00726             ipvol(2) = p2[2] + ipfac(0)*(p0[2]-p2[2]) + ipfac(1)*(p1[2]-p2[2]);
00727             break;
00728           }
00729         case ET_QUAD:
00730           {
00731             const POINT3D & p0 = points[faces[fnr][0]];
00732             const POINT3D & p1 = points[faces[fnr][1]];
00733             const POINT3D & p2 = points[faces[fnr][3]];
00734             ipvol(0) = p0[0] + ipfac(0)*(p1[0]-p0[0]) + ipfac(1)*(p2[0]-p0[0]);
00735             ipvol(1) = p0[1] + ipfac(0)*(p1[1]-p0[1]) + ipfac(1)*(p2[1]-p0[1]);
00736             ipvol(2) = p0[2] + ipfac(0)*(p1[2]-p0[2]) + ipfac(1)*(p2[2]-p0[2]);
00737             break;
00738           }
00739         default:
00740           throw Exception ("undefined facet type in Facet2ElementTrafo()\n");
00741 
00742         } 
00743       /*      cerr << "*** mapping integrationpoint for element " << eltype << " and facel " << fnr << " of type " << facettype << endl;
00744               cerr << "  * ipfac = " << ipfac;
00745               cerr << "  * ipvol = " << ipvol;*/
00746     }
00747     const IntegrationPoint operator()(const IntegrationPoint &ip1d) const 
00748     {
00749       IntegrationPoint elpoint;
00750       operator()(ip1d, elpoint);
00751       return elpoint;
00752     }
00753   };
00754 
00755 
00756 
00757 
00758 
00759 
00760 
00761 
00762 
00763 
00764 
00765 
00766 
00767   class NGS_DLL_HEADER BaseMappedIntegrationRule 
00768   { 
00769   protected:
00770     const IntegrationRule & ir;
00771     const ElementTransformation & eltrans;
00772     char * baseip;
00773     int incr;
00774 
00775   public:
00776     BaseMappedIntegrationRule (const IntegrationRule & air,
00777                                const ElementTransformation & aeltrans)
00778       : ir(air), eltrans(aeltrans) { ; }
00779 
00780     int Size() const { return ir.Size(); }
00781     const IntegrationRule & IR() const { return ir; }
00782     const ElementTransformation & GetTransformation () const { return eltrans; }
00783 
00784     const BaseMappedIntegrationPoint & operator[] (int i) const
00785     { return *static_cast<const BaseMappedIntegrationPoint*> ((void*)(baseip+i*incr)); }
00786   };
00787 
00788   template <int DIM_ELEMENT, int DIM_SPACE>
00789   class NGS_DLL_HEADER MappedIntegrationRule : public BaseMappedIntegrationRule
00790   {
00791     FlatArray< MappedIntegrationPoint<DIM_ELEMENT, DIM_SPACE> > sips;
00792   public:
00793     MappedIntegrationRule (const IntegrationRule & ir, 
00794                            const ElementTransformation & aeltrans, 
00795                            LocalHeap & lh);
00796 
00797     MappedIntegrationRule (const IntegrationRule & ir, 
00798                            const ElementTransformation & eltrans, 
00799                            int dummy,
00800                            LocalHeap & lh)
00801       : BaseMappedIntegrationRule (ir, eltrans), sips(ir.GetNIP(), lh)
00802     {
00803       baseip = (char*)(void*)(BaseMappedIntegrationPoint*)(&sips[0]);
00804       incr = (char*)(void*)(&sips[1]) - (char*)(void*)(&sips[0]);
00805     }
00806     
00807     MappedIntegrationPoint<DIM_ELEMENT, DIM_SPACE> & operator[] (int i) const
00808     { 
00809       return sips[i]; 
00810     }
00811   };
00812 
00813 #define SpecificIntegrationPoint MappedIntegrationPoint 
00814   /*
00815 #define BaseSpecificIntegrationPoint BaseMappedIntegrationPoint 
00816 #define DimSpecificIntegrationPoint DimMappedIntegrationPoint 
00817   */
00818 
00819 
00820 
00821 }
00822 
00823 
00824 #endif