/[escript]/branches/domexper/dudley/src/Quadrature.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Quadrature.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3169 - (show annotations)
Thu Sep 9 03:21:35 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 6056 byte(s)
Some moving

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Dudley: quadrature schemes */
18
19 /**************************************************************/
20
21 #include "Quadrature.h"
22
23
24 #define QUADNODES(_K_,_I_) quadNodes[INDEX2(_K_,_I_,DIM)]
25 #define QUADWEIGHTS(_I_) quadWeights[_I_]
26
27 /**************************************************************/
28
29 Dudley_QuadInfo Dudley_QuadInfoList[]={
30 {PointQuad, "Point", 0, 1, Dudley_Quad_getNodesPoint, Dudley_Quad_getNumNodesPoint} ,
31 {LineQuad, "Line", 1, 2, Dudley_Quad_getNodesLine, Dudley_Quad_getNumNodesLine} ,
32 {TriQuad, "Tri", 2, 3, Dudley_Quad_getNodesTri, Dudley_Quad_getNumNodesTri},
33 {TetQuad, "Tet", 3, 4, Dudley_Quad_getNodesTet, Dudley_Quad_getNumNodesTet},
34 {NoQuad, "NoType", 0, 1, Dudley_Quad_getNodesPoint, Dudley_Quad_getNumNodesPoint}
35 };
36
37 Dudley_QuadInfo* Dudley_QuadInfo_getInfo(Dudley_QuadTypeId id)
38 {
39 int ptr=0;
40 Dudley_QuadInfo* out=NULL;
41 while (Dudley_QuadInfoList[ptr].TypeId!=NoQuad && out==NULL) {
42 if (Dudley_QuadInfoList[ptr].TypeId==id) out=&(Dudley_QuadInfoList[ptr]);
43 ptr++;
44 }
45 if (out==NULL) {
46 Dudley_setError(VALUE_ERROR,"Dudley_QuadInfo_getInfo: canot find requested quadrature scheme.");
47 }
48 return out;
49 }
50
51 /**************************************************************/
52
53 /* get a quadrature scheme with numQuadNodes quadrature nodes for the tri */
54 /* as a queezed scheme on a quad [0,1]^2 */
55
56 void Dudley_Quad_getNodesTri(int numQuadNodes,double* quadNodes,double* quadWeights) {
57 // int i;
58 // double Q1,Q2,a,b,c,d,e,f,g,u,v,w;
59 #define DIM 2
60
61 /* the easy cases: */
62 if (numQuadNodes==1) {
63 QUADNODES(0,0)=1./3.;
64 QUADNODES(1,0)=1./3.;
65 QUADWEIGHTS(0)=1./2.;
66 } else if (numQuadNodes==3){
67 QUADNODES(0,0)=1./2.;
68 QUADNODES(1,0)=0.;
69 QUADWEIGHTS(0)=1./6.;
70 QUADNODES(0,1)=0.;
71 QUADNODES(1,1)=1./2.;
72 QUADWEIGHTS(1)=1./6.;
73 QUADNODES(0,2)=1./2.;
74 QUADNODES(1,2)=1./2.;
75 QUADWEIGHTS(2)=1./6.;
76 }
77 #undef DIM
78
79
80 }
81
82 /**************************************************************/
83
84 /* get a quadrature scheme with numQuadNodes quadrature nodes for the tet */
85 /* as a queezed scheme on a hex [0,1]^3 */
86
87 void Dudley_Quad_getNodesTet(int numQuadNodes,double* quadNodes,double* quadWeights) {
88 double alpha=0.58541019662496852;
89 double beta =0.1381966011250105;
90 #define DIM 3
91
92 /* the easy cases: */
93 if (numQuadNodes==1) {
94 QUADNODES(0,0)=0.25;
95 QUADNODES(1,0)=0.25;
96 QUADNODES(2,0)=0.25;
97 QUADWEIGHTS(0)=1./6.;
98 } else if (numQuadNodes==4){
99 QUADNODES(0,0)=beta;
100 QUADNODES(1,0)=beta;
101 QUADNODES(2,0)=beta;
102 QUADWEIGHTS(0)=1./24.;
103 QUADNODES(0,1)=alpha;
104 QUADNODES(1,1)=beta;
105 QUADNODES(2,1)=beta;
106 QUADWEIGHTS(1)=1./24.;
107 QUADNODES(0,2)=beta;
108 QUADNODES(1,2)=alpha;
109 QUADNODES(2,2)=beta;
110 QUADWEIGHTS(2)=1./24.;
111 QUADNODES(0,3)=beta;
112 QUADNODES(1,3)=beta;
113 QUADNODES(2,3)=alpha;
114 QUADWEIGHTS(3)=1./24.;
115 }
116 #undef DIM
117 }
118
119 /**************************************************************/
120
121 /* get a quadrature scheme with DSDV(2,2,i)= 0.;
122 DSDV(2,3,i)= 0.;
123 DSDV(3,1,i)= 0.;numQuadNodes quadrature nodes for a point. As there */
124 /* in no quadrature scheme for a point any value for numQuadNodes other than 0 throws */
125 /* an error. */
126
127 void Dudley_Quad_getNodesPoint(int numQuadNodes,double* quadNodes,double* quadWeights) {
128 if (numQuadNodes==0) {
129 return;
130 } else {
131 Dudley_setError(VALUE_ERROR,"Dudley_Quad_getNodesPoint: Illegal number of quadrature nodes.");
132 }
133 }
134
135 /**************************************************************/
136
137 /* get a quadrature scheme with numQuadNodes quadrature nodes on the line [0,1]: */
138 /* The nodes and weights are set from a table. */
139
140 void Dudley_Quad_getNodesLine(int numQuadNodes,double* quadNodes,double* quadWeights) {
141 switch(numQuadNodes) {
142 case 1:
143 quadNodes[0]=0.5;
144 quadWeights[0]=1.;
145 break;
146
147 case 2:
148 quadNodes[0]=(1.-.577350269189626)/2.;
149 quadNodes[1]=(1.+.577350269189626)/2.;
150 quadWeights[0]=.5;
151 quadWeights[1]=.5;
152 break;
153 default:
154 Dudley_setError(VALUE_ERROR,"Dudley_Quad_getNodesLine: Negative intergration order.");
155 break;
156 }
157 }
158
159
160 /**************************************************************/
161
162 /* The following functions Dudley_Quad_getNumNodes* return the nmber of quadrature points needed to */
163 /* achieve a certain accuracy. Notice that for Tet and Tri the the order is increased */
164 /* to consider the accuracy reduction through the construction process. */
165
166
167 int Dudley_Quad_getNumNodesPoint(int order) {
168 return 0;
169 }
170
171 int Dudley_Quad_getNumNodesLine(int order) {
172 char error_msg[LenErrorMsg_MAX];
173 if (order <0 ) {
174 Dudley_setError(VALUE_ERROR,"Dudley_Quad_getNumNodesPoint: Negative intergration order.");
175 return -1;
176 } else {
177 if (order > 2*MAX_numQuadNodesLine-1) {
178 sprintf(error_msg,"Dudley_Quad_getNumNodesPoint: requested integration order %d on line is too large (>%d).",
179 order,2*MAX_numQuadNodesLine-1);
180 Dudley_setError(VALUE_ERROR,error_msg);
181 return -1;
182 } else {
183 Dudley_resetError();
184 return order/2+1;
185 }
186 }
187 }
188
189 int Dudley_Quad_getNumNodesTri(int order) {
190 if (order<=1) {
191 return 1;
192 } else if (order<=2){
193 return 3;
194 } else {
195 return -1;
196 }
197 }
198
199 int Dudley_Quad_getNumNodesTet(int order) {
200 if (order<=1) {
201 return 1;
202 } else if (order<=2){
203 return 4;
204 } else {
205 return -1;
206 }
207 }
208

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26