← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Geo/GaussQuadratureHex.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2009-02-17 10:12:27 UTC
  • mfrom: (1.2.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20090217101227-mdrolkldak2pgd2i
Tags: 2.3.0.dfsg-1
* New upstream release
  + major graphics and GUI code refactoring; 
  + new full-quad/hexa subdivision algorithm (removed 
    Mesh.RecombineAlgo);
  + improved automatic transfinite corner selection (now also 
    for volumes); 
  + improved visibility browser; new automatic adaptive visualization
    for high-order simplices;
  + modified arrow size, clipping planes and transform options; many
    improvements and
  + bug fixes all over the place.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Geo/GaussQuadratureHex.cpp
1
 
// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
2
 
//
3
 
// See the LICENSE.txt file for license information. Please report all
4
 
// bugs and problems to <gmsh@geuz.org>.
5
 
 
6
 
#include "MElement.h"
7
 
#include "GaussLegendreSimplex.h"
8
 
#include "GaussLegendre1D.h"
9
 
 
10
 
const double a1  = 0.40824826;
11
 
const double ma1 =-0.40824826;
12
 
const double a2  = 0.81649658;
13
 
const double ma2 =-0.81649658;
14
 
const double b1  = 0.70710678;
15
 
const double mb1 =-0.70710678;
16
 
const double c1  = 0.57735027;
17
 
const double mc1 =-0.57735027;
18
 
const double w1  = 1.3333333333;
19
 
const double xh6[6] = { a1,  a1, ma1, ma1, ma2,  a2};
20
 
const double yh6[6] = { b1, mb1,  b1, mb1,  0.,  0.};
21
 
const double zh6[6] = {mc1, mc1,  c1,  c1, mc1,  c1};
22
 
const double ph6[6] = { w1,  w1,  w1,  w1,  w1,  w1};
23
 
 
24
 
IntPt GQH1[1] = 
25
 
{
26
 
  { {0.0,0.0,0.0},8.0 }
27
 
};
28
 
 
29
 
IntPt GQH6[6] = 
30
 
{
31
 
  { {xh6[0],yh6[0],zh6[0]},ph6[0] },
32
 
  { {xh6[1],yh6[1],zh6[1]},ph6[1] },
33
 
  { {xh6[2],yh6[2],zh6[2]},ph6[2] },
34
 
  { {xh6[3],yh6[3],zh6[3]},ph6[3] },
35
 
  { {xh6[4],yh6[4],zh6[4]},ph6[4] },
36
 
  { {xh6[5],yh6[5],zh6[5]},ph6[5] }
37
 
};
38
 
 
39
 
const double xh8[8] = {0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626,
40
 
                   0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626};
41
 
const double yh8[8] = {0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626,
42
 
                   0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626};
43
 
 
44
 
const double zh8[8] = {-0.577350269189626,-0.577350269189626,-0.577350269189626,-0.577350269189626,
45
 
                   0.577350269189626,0.577350269189626,0.577350269189626,0.577350269189626};
46
 
const double ph8[8] = {1.,1.,1.,1.,1.,1.,1.,1.};
47
 
 
48
 
IntPt GQH8[8] = 
49
 
{
50
 
  { {xh8[0],yh8[0],zh8[0]},ph8[0] },
51
 
  { {xh8[1],yh8[1],zh8[1]},ph8[1] },
52
 
  { {xh8[2],yh8[2],zh8[2]},ph8[2] },
53
 
  { {xh8[3],yh8[3],zh8[3]},ph8[3] },
54
 
  { {xh8[4],yh8[4],zh8[4]},ph8[4] },
55
 
  { {xh8[5],yh8[5],zh8[5]},ph8[5] },
56
 
  { {xh8[6],yh8[6],zh8[6]},ph8[6] },
57
 
  { {xh8[7],yh8[7],zh8[7]},ph8[7] }
58
 
};
59
 
 
60
 
 
61
 
IntPt *getGQHPts(int order);
62
 
int getNGQHPts(int order);
63
 
 
64
 
IntPt * GQH[17] = {GQH1,GQH1,GQH6,GQH8,0,0,0,0,0,0,0,0,0,0,0,0,0};
65
 
int GQHnPt[4] = {1,4,5,15};
66
 
 
67
 
IntPt *getGQHPts(int order)
68
 
{ 
69
 
 
70
 
  if(order<2)return GQH[order];
71
 
  if(order == 2)return GQH[3]; 
72
 
  if(order == 3)return GQH[3]; 
73
 
  
74
 
  int n = (order+3)/2;
75
 
  int index = n-2 + 4;
76
 
  if(!GQH[index])
77
 
    {
78
 
      double *pt,*wt;
79
 
      //      printf("o = %d n = %d i = %d\n",order,n*n*n,index);
80
 
      gmshGaussLegendre1D(n,&pt,&wt);
81
 
      GQH[index] = new IntPt[n*n*n];
82
 
      int l = 0;
83
 
      for(int i=0; i < n; i++) {
84
 
        for(int j=0; j < n; j++) {
85
 
          for(int k=0; k < n; k++) {
86
 
            GQH[index][l].pt[0] = pt[i];
87
 
            GQH[index][l].pt[1] = pt[j];
88
 
            GQH[index][l].pt[2] = pt[k];
89
 
            GQH[index][l++].weight = wt[i]*wt[j]*wt[k];
90
 
            //      printf ("%f %f %f %f\n",pt[i],pt[j],pt[k],wt[i]*wt[j]*wt[k]);
91
 
          }
92
 
        }
93
 
      }
94
 
    }
95
 
  return GQH[index];
96
 
}
97
 
 
98
 
int getNGQHPts(int order)
99
 
{ 
100
 
  if(order == 3)return 8;
101
 
  if(order == 2)return 8;
102
 
  if(order < 2)
103
 
    return GQHnPt[order]; 
104
 
  return ((order+3)/2)*((order+3)/2)*((order+3)/2);
105
 
}
106
 
 
107
 
 
108
 
 
109
 
 
110