NGSolve
4.9
|
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