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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3114 - (show annotations)
Fri Aug 27 05:26:25 2010 UTC (9 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 18375 byte(s)
It doesn't pass all tests but this is major progress

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: Shape functions */
18
19 /**************************************************************/
20
21 #include "ShapeFunctions.h"
22
23
24
25 Dudley_ShapeFunctionInfo Dudley_ShapeFunction_InfoList[]={
26 {Point1Shape, "Point1", 0, 1, 1, 1, Dudley_Shape_Point1 } ,
27 {Line2Shape, "Line2", 1, 2, 1, 2, Dudley_Shape_Line2 } ,
28 {Line3Shape, "Line3", 1, 3, 2, 2, Dudley_Shape_Line3 },
29 {Line4Shape, "Line4", 1, 4, 3, 2, Dudley_Shape_Line4 },
30 {Tri3Shape, "Tri3", 2, 3, 1, 3, Dudley_Shape_Tri3 },
31 {Tri6Shape, "Tri6", 2, 6, 2, 3, Dudley_Shape_Tri6 },
32 {Tri9Shape, "Tri9", 2, 9, 3, 3, Dudley_Shape_Tri9 },
33 {Tri10Shape, "Tri10", 2, 10, 3, 3, Dudley_Shape_Tri10, },
34 {Tet4Shape, "Tet4", 3, 4, 1, 4, Dudley_Shape_Tet4, },
35 {Tet10Shape, "Tet10", 3, 10, 2, 4, Dudley_Shape_Tet10, },
36 {Tet16Shape, "Tet16", 3, 16, 3, 4, Dudley_Shape_Tet16, },
37 {NoShape, "NoType", 0, 1, 1, 1, Dudley_Shape_Point1 }
38 };
39
40
41 /******************************************************************************************************************************
42
43 creates an evaluation of the ShapeFunction on the given quadrature scheme.
44 if the spatial dimension of the scheme and the shape functions don't match
45
46 if QuadNodes==Null or QuadWeights==Null the shape functions method is used to generate a quadrature scheme with numQuasNodes
47 nodes. otherwise its assumed that a quadraure scheme is given on these array and copy is created within the structure.
48
49 */
50 Dudley_ShapeFunction* Dudley_ShapeFunction_alloc(Dudley_ShapeFunctionTypeId id,int numQuadDim, int numQuadNodes, double *QuadNodes, double *QuadWeights) {
51 Dudley_ShapeFunction *out=NULL;
52 int numDim, numShapes, i, q;
53
54 numDim=Dudley_ShapeFunction_InfoList[id].numDim;
55 numShapes=Dudley_ShapeFunction_InfoList[id].numShapes;
56
57 if (numQuadDim>numDim) {
58 Dudley_setError(VALUE_ERROR,"Dudley_ShapeFunction_alloc: spatial dimension of quadrature scheme is bigger then spatial dimension of shape function.");
59 return NULL;
60 }
61
62 /* allocate the Dudley_ShapeFunction to be returned: */
63
64 out=MEMALLOC(1,Dudley_ShapeFunction);
65 if (Dudley_checkPtr(out)) return NULL;
66
67
68 out->Type=Dudley_ShapeFunction_getInfo(id);
69 out->numQuadNodes=numQuadNodes;
70 out->QuadNodes=NULL;
71 out->QuadWeights=NULL;
72 out->S=NULL;
73 out->dSdv=NULL;
74 out->reference_counter=0;
75
76 /* allocate memory: */
77
78 out->QuadNodes=MEMALLOC(numQuadNodes*numDim,double);
79 out->QuadWeights=MEMALLOC(numQuadNodes,double);
80 out->S=MEMALLOC(numShapes*numQuadNodes,double);
81 out->dSdv=MEMALLOC(numShapes*numDim*numQuadNodes,double);
82 if ( Dudley_checkPtr(out->QuadNodes) || Dudley_checkPtr(out->QuadWeights) || Dudley_checkPtr(out->S) || Dudley_checkPtr(out->dSdv) ) {
83 Dudley_ShapeFunction_dealloc(out);
84 return NULL;
85 }
86
87 /* set the quadrature nodes (missing values are filled with 0): */
88
89 for (q=0;q<numQuadNodes;q++) {
90 for (i=0;i<numQuadDim;i++) out->QuadNodes[INDEX2(i,q,numDim)]=QuadNodes[INDEX2(i,q,numQuadDim)];
91 for (i=numQuadDim;i<numDim;i++) out->QuadNodes[INDEX2(i,q,numDim)]=0;
92 out->QuadWeights[q]=QuadWeights[q];
93 }
94
95 /* eval shape functions on quadrature node: */
96
97 out->Type->getValues(numQuadNodes,out->QuadNodes,out->S,out->dSdv);
98
99 if (! Dudley_noError()) {
100 Dudley_ShapeFunction_dealloc(out);
101 return NULL;
102 }
103
104 /* all done: */
105 out->reference_counter=1;
106 return out;
107 }
108
109 Dudley_ShapeFunction* Dudley_ShapeFunction_reference(Dudley_ShapeFunction* in) {
110 if (in!=NULL) ++(in->reference_counter);
111 return in;
112 }
113 /**************************************************************/
114
115 void Dudley_ShapeFunction_dealloc(Dudley_ShapeFunction* in) {
116 if (in!=NULL) {
117 in->reference_counter--;
118 if (in->reference_counter<1) {
119 MEMFREE(in->QuadNodes);
120 MEMFREE(in->QuadWeights);
121 MEMFREE(in->S);
122 MEMFREE(in->dSdv);
123 MEMFREE(in);
124 }
125 }
126 }
127
128 /**************************************************************/
129
130 Dudley_ShapeFunctionTypeId Dudley_ShapeFunction_getTypeId(char* element_type)
131 {
132 int ptr=0;
133 Dudley_ShapeFunctionTypeId out=NoShape;
134 while (Dudley_ShapeFunction_InfoList[ptr].TypeId!=NoShape && out==NoShape) {
135 if (strcmp(element_type,Dudley_ShapeFunction_InfoList[ptr].Name)==0) out=Dudley_ShapeFunction_InfoList[ptr].TypeId;
136 ptr++;
137 }
138 return out;
139 }
140
141 Dudley_ShapeFunctionInfo* Dudley_ShapeFunction_getInfo(Dudley_ShapeFunctionTypeId id)
142 {
143 int ptr=0;
144 Dudley_ShapeFunctionInfo* out=NULL;
145 while (Dudley_ShapeFunction_InfoList[ptr].TypeId!=NoShape && out==NULL) {
146 if (Dudley_ShapeFunction_InfoList[ptr].TypeId==id) out=&(Dudley_ShapeFunction_InfoList[ptr]);
147 ptr++;
148 }
149 if (out==NULL) {
150 Dudley_setError(VALUE_ERROR,"Dudley_ShapeFunctionInfo_getInfo: canot find requested shape function");
151 }
152 return out;
153 }
154
155 /**************************************************************/
156
157 #define V(_K_,_I_) v[INDEX2((_K_)-1,(_I_),DIM)]
158 #define S(_J_,_I_) s[S_INDEX((_J_)-1,(_I_),NUMSHAPES)]
159 #define DSDV(_J_,_K_,_I_) dsdv[DSDV_INDEX((_J_)-1,(_K_)-1,(_I_),NUMSHAPES,DIM)]
160
161 /**************************************************************/
162
163 void Dudley_Shape_Point1(int NumV,double* v,double* s,double* dsdv) {
164 #define NUMSHAPES 1
165 #define DIM 0
166 int i;
167 #pragma ivdep
168 for (i=0;i<NumV;i++) {
169 S(1,i)=1.;
170 }
171 #undef NUMSHAPES
172 #undef DIM
173 }
174
175 /**************************************************************/
176
177 void Dudley_Shape_Line2(int NumV,double* v,double* s,double* dsdv) {
178 #define NUMSHAPES 2
179 #define DIM 1
180 register double x;
181 int i;
182 #pragma ivdep
183 for (i=0;i<NumV;i++) {
184 x=V(1,i);
185 S(1,i)=1.-x;
186 S(2,i)= x;
187 DSDV(1,1,i)=-1.;
188 DSDV(2,1,i)= 1.;
189 }
190 #undef NUMSHAPES
191 #undef DIM
192 }
193
194 /**************************************************************/
195
196 void Dudley_Shape_Line3(int NumV,double* v,double* s,double* dsdv) {
197 #define NUMSHAPES 3
198 #define DIM 1
199 register double x;
200 int i;
201 #pragma ivdep
202 for (i=0;i<NumV;i++) {
203 x=V(1,i);
204 S(1,i)=(2.*x -1. )*(x -1.);
205 S(2,i)=(2.*x -1.)*x;
206 S(3,i)= 4.*x*(1. -x );
207 DSDV(1,1,i)= 4.*x -3.;
208 DSDV(2,1,i)= 4.*x -1.;
209 DSDV(3,1,i)=-8.*x+4.;
210 }
211 #undef NUMSHAPES
212 #undef DIM
213 }
214
215 /**************************************************************/
216
217 void Dudley_Shape_Line4(int NumV,double* v,double* s,double* dsdv) {
218 #define NUMSHAPES 4
219 #define DIM 1
220 register double x;
221 int i;
222 #pragma ivdep
223 for (i=0;i<NumV;i++) {
224 x=V(1,i);
225 S(1,i)=(10.)+(-5.5)*x+(9.)*x*x+(-4.5)*x*x*x ;
226 S(2,i)=(10.)*x+(-4.5)*x*x+(4.5)*x*x*x ;
227 S(3,i)=(9.)*x+(-22.5)*x*x+(13.5)*x*x*x ;
228 S(4,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x;
229 DSDV(1,1,i)=(-5.5)+(18.)*x+(-13.5)*x*x;
230 DSDV(2,1,i)=(10.)+(-9.)*x+(13.5)*x*x;
231 DSDV(3,1,i)=(9.)+(-45.)*x+(0.405e2)*x*x;
232 DSDV(4,1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x;
233 }
234 #undef NUMSHAPES
235 #undef DIM
236 }
237
238 /**************************************************************/
239
240 void Dudley_Shape_Tri3(int NumV,double* v,double* s,double* dsdv) {
241 #define NUMSHAPES 3
242 #define DIM 2
243 register double x,y;
244 int i;
245 #pragma ivdep
246 for (i=0;i<NumV;i++) {
247 x=V(1,i);
248 y=V(2,i);
249 S(1,i)=1.-x-y;
250 S(2,i)= x;
251 S(3,i)= y;
252 DSDV(1,1,i)=-1.;
253 DSDV(1,2,i)=-1.;
254 DSDV(2,1,i)= 1.;
255 DSDV(2,2,i)= 0.;
256 DSDV(3,1,i)= 0.;
257 DSDV(3,2,i)= 1.;
258 }
259 #undef NUMSHAPES
260 #undef DIM
261 }
262
263 /**************************************************************/
264
265 void Dudley_Shape_Tri6(int NumV,double* v,double* s,double* dsdv) {
266 #define NUMSHAPES 6
267 #define DIM 2
268 register double x,y;
269 int i;
270 #pragma ivdep
271 for (i=0;i<NumV;i++) {
272 x=V(1,i);
273 y=V(2,i);
274 S(1,i)= (1. -x -y)*(1. -2.*x -2.* y);
275 S(2,i)= x*(2.* x -1.);
276 S(3,i)= y*(2.* y -1.);
277 S(4,i)= (1. -x -y)*4.* x;
278 S(5,i)= 4.*x*y;
279 S(6,i)= (1. -x -y)*4.* y;
280 DSDV(1,1,i)= -3.+4.*x+4.*y;
281 DSDV(1,2,i)= -3.+4.*x+4.*y;
282 DSDV(2,1,i)= -1.+4.*x;
283 DSDV(2,2,i)= 0.;
284 DSDV(3,1,i)= 0.;
285 DSDV(3,2,i)= -1. +4.*y;
286 DSDV(4,1,i)= 4. -8.*x -4.*y;
287 DSDV(4,2,i)= -4.*x;
288 DSDV(5,1,i)= 4.*y;
289 DSDV(5,2,i)= 4.*x;
290 DSDV(6,1,i)= -4.*y;
291 DSDV(6,2,i)= 4. -4.*x -8.*y;
292 }
293 #undef NUMSHAPES
294 #undef DIM
295 }
296
297 /**************************************************************/
298
299 void Dudley_Shape_Tri9(int NumV,double* v,double* s,double* dsdv) {
300 #define NUMSHAPES 9
301 #define DIM 2
302 register double x,y;
303 int i;
304 #pragma ivdep
305 for (i=0;i<NumV;i++) {
306 x=V(1,i);
307 y=V(2,i);
308 S(1,i)=(10.)+(-5.5)*x+(-5.5)*y+(9.)*x*x+(-4.5)*x*x*x+(9.)*y*y+(-4.5)*y*y*y+(4.5)*x*y*y+(4.5)*x*x*y;
309 S(2,i)=(10.)*x+(-4.5)*x*x+(4.5)*x*x*x;
310 S(3,i)=(10.)*y+(-4.5)*y*y+(4.5)*y*y*y;
311 S(4,i)=(9.)*x+(-22.5)*x*x+(13.5)*x*x*x+(-9.)*x*y*y+(4.5)*x*x*y;
312 S(5,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x+(4.5)*x*y*y+(-9.)*x*x*y;
313 S(6,i)=(-4.5)*x*y*y+(9.)*x*x*y;
314 S(7,i)=(9.)*x*y*y+(-4.5)*x*x*y;
315 S(8,i)=(-4.5)*y+(18.)*y*y+(-13.5)*y*y*y+(-9.)*x*y*y+(4.5)*x*x*y;
316 S(9,i)=(9.)*y+(-22.5)*y*y+(13.5)*y*y*y+(4.5)*x*y*y+(-9.)*x*x*y;
317 DSDV(1, 1,i)=(-5.5)+(18.)*x+(-13.5)*x*x+(4.5)*y*y+(9.)*x*y;
318 DSDV(2, 1,i)=(10.)+(-9.)*x+(13.5)*x*x;
319 DSDV(3, 1,i)= 0.;
320 DSDV(4, 1,i)=(9.)+(-45.)*x+(0.405e2)*x*x+(-9.)*y*y+(9.)*x*y;
321 DSDV(5, 1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x+(4.5)*y*y+(-18.)*x*y;
322 DSDV(6, 1,i)=(-4.5)*y*y+(18.)*x*y;
323 DSDV(7, 1,i)=(9.)*y*y+(-9.)*x*y;
324 DSDV(8, 1,i)=(-9.)*y*y+(9.)*x*y;
325 DSDV(9, 1,i)=(4.5)*y*y+(-18.)*x*y;
326 DSDV(1, 2,i)=(-5.5)+(18.)*y+(-13.5)*y*y+(9.)*x*y+(4.5)*x*x;
327 DSDV(2, 2,i)= 0.;
328 DSDV(3, 2,i)=(10.)+(-9.)*y+(13.5)*y*y;
329 DSDV(4, 2,i)=(-18.)*x*y+(4.5)*x*x;
330 DSDV(5, 2,i)=(9.)*x*y+(-9.)*x*x;
331 DSDV(6, 2,i)=(-9.)*x*y+(9.)*x*x;
332 DSDV(7, 2,i)=(18.)*x*y+(-4.5)*x*x;
333 DSDV(8, 2,i)=(-4.5)+(36.)*y+(-0.405e2)*y*y+(-18.)*x*y+(4.5)*x*x;
334 DSDV(9, 2,i)=(9.)+(-45.)*y+(0.405e2)*y*y+(9.)*x*y+(-9.)*x*x;
335 }
336 #undef NUMSHAPES
337 #undef DIM
338 }
339
340 /**************************************************************/
341
342 void Dudley_Shape_Tri10(int NumV,double* v,double* s,double* dsdv) {
343 #define NUMSHAPES 10
344 #define DIM 2
345 register double x,y;
346 int i;
347 #pragma ivdep
348 for (i=0;i<NumV;i++) {
349 x=V(1,i);
350 y=V(2,i);
351 S(1,i)=(10.)+(-5.5)*x+(-5.5)*y+(9.)*x*x+(-4.5)*x*x*x+(9.)*y*y+(-4.5)*y*y*y+(-13.5)*x*y*y+(-13.5)*x*x*y+(18.)*x*y;
352 S(2,i)=(10.)*x+(-4.5)*x*x+(4.5)*x*x*x;
353 S(3,i)=(10.)*y+(-4.5)*y*y+(4.5)*y*y*y;
354 S(4,i)=(9.)*x+(-22.5)*x*x+(13.5)*x*x*x+(13.5)*x*y*y+(0.27e2)*x*x*y+(-22.5)*x*y;
355 S(5,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x+(-13.5)*x*x*y+(4.5)*x*y;
356 S(6,i)=(13.5)*x*x*y+(-4.5)*x*y;
357 S(7,i)=(13.5)*x*y*y+(-4.5)*x*y;
358 S(8,i)=(-4.5)*y+(18.)*y*y+(-13.5)*y*y*y+(-13.5)*x*y*y+(4.5)*x*y;
359 S(9,i)=(9.)*y+(-22.5)*y*y+(13.5)*y*y*y+(0.27e2)*x*y*y+(13.5)*x*x*y+(-22.5)*x*y;
360 S(10,i)=(-0.27e2)*x*y*y+(-0.27e2)*x*x*y+(0.27e2)*x*y;
361 DSDV(1, 1,i)=(-5.5)+(18.)*x+(-13.5)*x*x+(-13.5)*y*y+(-0.27e2)*x*y+(18.)*y;
362 DSDV(2, 1,i)=(10.)+(-9.)*x+(13.5)*x*x;
363 DSDV(3, 1,i)= 0.;
364 DSDV(4, 1,i)=(9.)+(-45.)*x+(0.405e2)*x*x+(13.5)*y*y+(0.54e2)*x*y+(-22.5)*y;
365 DSDV(5, 1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x+(-0.27e2)*x*y+(4.5)*y;
366 DSDV(6, 1,i)=(0.27e2)*x*y+(-4.5)*y;
367 DSDV(7, 1,i)=(13.5)*y*y+(-4.5)*y;
368 DSDV(8, 1,i)=(-13.5)*y*y+(4.5)*y;
369 DSDV(9, 1,i)=(0.27e2)*y*y+(0.27e2)*x*y+(-22.5)*y;
370 DSDV(10, 1,i)=(-0.27e2)*y*y+(-0.54e2)*x*y+(0.27e2)*y;
371 DSDV(1, 2,i)=(-5.5)+(18.)*y+(-13.5)*y*y+(-0.27e2)*x*y+(-13.5)*x*x+(18.)*x;
372 DSDV(2, 2,i)=0.;
373 DSDV(3, 2,i)=(10.)+(-9.)*y+(13.5)*y*y;
374 DSDV(4, 2,i)=(0.27e2)*x*y+(0.27e2)*x*x+(-22.5)*x;
375 DSDV(5, 2,i)=(-13.5)*x*x+(4.5)*x;
376 DSDV(6, 2,i)=(13.5)*x*x+(-4.5)*x;
377 DSDV(7, 2,i)=(0.27e2)*x*y+(-4.5)*x;
378 DSDV(8, 2,i)=(-4.5)+(36.)*y+(-0.405e2)*y*y+(-0.27e2)*x*y+(4.5)*x;
379 DSDV(9, 2,i)=(9.)+(-45.)*y+(0.405e2)*y*y+(0.54e2)*x*y+(13.5)*x*x+(-22.5)*x;
380 DSDV(10, 2,i)=(-0.54e2)*x*y+(-0.27e2)*x*x+(0.27e2)*x;
381 }
382 #undef NUMSHAPES
383 #undef DIM
384 }
385
386
387 /**************************************************************/
388
389 void Dudley_Shape_Tet4(int NumV,double* v,double* s,double* dsdv) {
390 #define NUMSHAPES 4
391 #define DIM 3
392 register double x,y,z;
393 int i;
394 #pragma ivdep
395 for (i=0;i<NumV;i++) {
396 x=V(1,i);
397 y=V(2,i);
398 z=V(3,i);
399 S(1,i)=1.-x-y-z;
400 S(2,i)=x;
401 S(3,i)=y;
402 S(4,i)=z;
403 DSDV(1,1,i)=-1.;
404 DSDV(1,2,i)=-1.;
405 DSDV(1,3,i)=-1.;
406 DSDV(2,1,i)= 1.;
407 DSDV(2,2,i)= 0.;
408 DSDV(2,3,i)= 0.;
409 DSDV(3,1,i)= 0.;
410 DSDV(3,2,i)= 1.;
411 DSDV(3,3,i)= 0.;
412 DSDV(4,1,i)= 0.;
413 DSDV(4,2,i)= 0.;
414 DSDV(4,3,i)= 1.;
415 }
416 #undef NUMSHAPES
417 #undef DIM
418 }
419
420 /**************************************************************/
421
422 void Dudley_Shape_Tet10(int NumV,double* v,double* s,double* dsdv) {
423 #define NUMSHAPES 10
424 #define DIM 3
425 register double x,y,z;
426 int i;
427 #pragma ivdep
428 for (i=0;i<NumV;i++) {
429 x=V(1,i);
430 y=V(2,i);
431 z=V(3,i);
432 S(1,i) = (1.-x-y-z)*(1.-2.*x-2.*y-2.*z);
433 S(2,i) = x*(2.*x-1.);
434 S(3,i) = y*(2.*y-1.);
435 S(4,i) = z*(2.*z-1.);
436 S(5,i) = (1.-x-y-z)*4.*x;
437 S(6,i) = 4.*x*y;
438 S(7,i) = (1.-x-y-z)*4.*y;
439 S(8,i) = (1.-x-y-z)*4.*z;
440 S(9,i) = 4.*x*z;
441 S(10,i)= 4.*y*z;
442
443 DSDV(1,1,i)= -3.+4.*x+4.*y+4.*z;
444 DSDV(1,2,i)= -3.+4.*x+4.*y+4.*z;
445 DSDV(1,3,i)= -3.+4.*x+4.*y+4.*z;
446
447
448 DSDV(2,1,i)= -1.+4.*x;
449 DSDV(2,2,i)= 0.;
450 DSDV(2,3,i)= 0.;
451
452 DSDV(3,1,i)= 0.;
453 DSDV(3,2,i)= -1. +4.*y;
454 DSDV(3,3,i)= 0.;
455
456 DSDV(4,1,i)= 0.;
457 DSDV(4,2,i)= 0.;
458 DSDV(4,3,i)= -1. +4.*z;
459
460 DSDV(5,1,i)= 4. -8.*x -4.*y -4.*z;
461 DSDV(5,2,i)= -4.*x;
462 DSDV(5,3,i)= -4.*x;
463
464 DSDV(6,1,i)= 4.*y;
465 DSDV(6,2,i)= 4.*x;
466 DSDV(6,3,i)= 0.;
467
468 DSDV(7,1,i)= -4.*y;
469 DSDV(7,2,i)= 4. -4.*x -8.*y -4.*z;
470 DSDV(7,3,i)= -4.*y;
471
472 DSDV(8,1,i)= -4.*z;
473 DSDV(8,2,i)= -4.*z;
474 DSDV(8,3,i)= 4. -4.*x -4.*y -8.*z;
475
476 DSDV(9,1,i)= 4.*z;
477 DSDV(9,2,i)= 0.;
478 DSDV(9,3,i)= 4.*x;
479
480 DSDV(10,1,i)= 0.;
481 DSDV(10,2,i)= 4.*z;
482 DSDV(10,3,i)= 4.*y;
483 }
484 #undef NUMSHAPES
485 #undef DIM
486 }
487
488 /**************************************************************/
489
490 void Dudley_Shape_Tet16(int NumV,double* v,double* s,double* dsdv) {
491 #define NUMSHAPES 16
492 #define DIM 3
493 register double x,y,z;
494 int i;
495 #pragma ivdep
496 for (i=0;i<NumV;i++) {
497 x=V(1,i);
498 y=V(2,i);
499 z=V(3,i);
500 S(1,i)=(10.)+(-5.5)*x+(-5.5)*y+(-5.5)*z+(9.)*x*x+(-4.5)*x*x*x+(4.5)*x*x*y+(4.5)*x*y*y+(-4.5)*y*y*y+(9.)*y*y+(9.)*z*z+(4.5)*x*x*z+(4.5)*y*y*z+(-4.5)*z*z*z+(4.5)*x*z*z+(4.5)*y*z*z;
501 S(2,i)=(1.e0)*x+(-4.5)*x*x+(4.5)*x*x*x;
502 S(3,i)=(1.e0)*y+(4.5)*y*y*y+(-4.5)*y*y;
503 S(4,i)=(1.e0)*z+(-4.5)*z*z+(4.5)*z*z*z;
504 S(5,i)=(9.)*x+(-22.5)*x*x+(13.5)*x*x*x+(4.5)*x*x*y+(-9.)*x*y*y+(4.5)*x*x*z+(-9.)*x*z*z;
505 S(6,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x+(-9.)*x*x*y+(4.5)*x*y*y+(-9.)*x*x*z+(4.5)*x*z*z;
506 S(7,i)=(9.)*x*x*y+(-4.5)*x*y*y;
507 S(8,i)=(-4.5)*x*x*y+(9.)*x*y*y;
508 S(9,i)=(-4.5)*y+(4.5)*x*x*y+(-9.)*x*y*y+(-13.5)*y*y*y+(18.)*y*y+(-9.)*y*y*z+(4.5)*y*z*z;
509 S(10,i)=(9.)*y+(-9.)*x*x*y+(4.5)*x*y*y+(13.5)*y*y*y+(-22.5)*y*y+(4.5)*y*y*z+(-9.)*y*z*z;
510 S(11,i)=(9.)*z+(-22.5)*z*z+(-9.)*x*x*z+(-9.)*y*y*z+(13.5)*z*z*z+(4.5)*x*z*z+(4.5)*y*z*z;
511 S(12,i)=(9.)*x*x*z+(-4.5)*x*z*z;
512 S(13,i)=(9.)*y*y*z+(-4.5)*y*z*z;
513 S(14,i)=(-4.5)*z+(18.)*z*z+(4.5)*x*x*z+(4.5)*y*y*z+(-13.5)*z*z*z+(-9.)*x*z*z+(-9.)*y*z*z;
514 S(15,i)=(-4.5)*x*x*z+(9.)*x*z*z;
515 S(16,i)=(-4.5)*y*y*z+(9.)*y*z*z;
516 DSDV(1, 1,i)=(-5.5)+(18.)*x+(-13.5)*x*x+(9.)*x*y+(4.5)*y*y+(9.)*x*z+(4.5)*z*z;
517 DSDV(2, 1,i)=(1.e0)+(-9.)*x+(13.5)*x*x;
518 DSDV(3, 1,i)= 0.;
519 DSDV(4, 1,i)= 0.;
520 DSDV(5, 1,i)=(9.)+(-45.)*x+(0.405e2)*x*x+(9.)*x*y+(-9.)*y*y+(9.)*x*z+(-9.)*z*z;
521 DSDV(6, 1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x+(-18.)*x*y+(4.5)*y*y+(-18.)*x*z+(4.5)*z*z;
522 DSDV(7, 1,i)=(18.)*x*y+(-4.5)*y*y;
523 DSDV(8, 1,i)=(-9.)*x*y+(9.)*y*y;
524 DSDV(9, 1,i)=(9.)*x*y+(-9.)*y*y;
525 DSDV(10, 1,i)=(-18.)*x*y+(4.5)*y*y;
526 DSDV(11, 1,i)=(-18.)*x*z+(4.5)*z*z;
527 DSDV(12, 1,i)=(18.)*x*z+(-4.5)*z*z;
528 DSDV(13, 1,i)=0.;
529 DSDV(14, 1,i)=(9.)*x*z+(-9.)*z*z;
530 DSDV(15, 1,i)=(-9.)*x*z+(9.)*z*z;
531 DSDV(16, 1,i)=0.;
532 DSDV(1, 2,i)=(-5.5)+(4.5)*x*x+(9.)*x*y+(-13.5)*y*y+(18.)*y+(9.)*y*z+(4.5)*z*z;
533 DSDV(2, 2,i)=0.;
534 DSDV(3, 2,i)=(1.e0)+(13.5)*y*y+(-9.)*y;
535 DSDV(4, 2,i)=0.;
536 DSDV(5, 2,i)=(4.5)*x*x+(-18.)*x*y;
537 DSDV(6, 2,i)=(-9.)*x*x+(9.)*x*y;
538 DSDV(7, 2,i)=(9.)*x*x+(-9.)*x*y;
539 DSDV(8, 2,i)=(-4.5)*x*x+(18.)*x*y;
540 DSDV(9, 2,i)=(-4.5)+(4.5)*x*x+(-18.)*x*y+(-0.405e2)*y*y+(36.)*y+(-18.)*y*z+(4.5)*z*z;
541 DSDV(10, 2,i)=(9.)+(-9.)*x*x+(9.)*x*y+(0.405e2)*y*y+(-45.)*y+(9.)*y*z+(-9.)*z*z;
542 DSDV(11, 2,i)=(-18.)*y*z+(4.5)*z*z;
543 DSDV(12, 2,i)=0.;
544 DSDV(13, 2,i)=(18.)*y*z+(-4.5)*z*z;
545 DSDV(14, 2,i)=(9.)*y*z+(-9.)*z*z;
546 DSDV(15, 2,i)=0.;
547 DSDV(16, 2,i)=(-9.)*y*z+(9.)*z*z;
548 DSDV(1, 3,i)=(-5.5)+(18.)*z+(4.5)*x*x+(4.5)*y*y+(-13.5)*z*z+(.9e1)*x*z+(9.)*y*z;
549 DSDV(2, 3,i)= 0.;
550 DSDV(3, 3,i)= 0.;
551 DSDV(4, 3,i)=(1.e0)+(-9.)*z+(13.5)*z*z;
552 DSDV(5, 3,i)=(4.5)*x*x+(-18.)*x*z;
553 DSDV(6, 3,i)=(-9.)*x*x+(9.)*x*z;
554 DSDV(7, 3,i)= 0.;
555 DSDV(8, 3,i)= 0.;
556 DSDV(9, 3,i)=(-9.)*y*y+(9.)*y*z;
557 DSDV(10, 3,i)=(4.5)*y*y+(-18.)*y*z;
558 DSDV(11, 3,i)=(9.)+(-45.)*z+(-9.)*x*x+(-9.)*y*y+(0.405e2)*z*z+(.9e1)*x*z+(9.)*y*z;
559 DSDV(12, 3,i)=(9.)*x*x+(-9.)*x*z;
560 DSDV(13, 3,i)=(9.)*y*y+(-9.)*y*z;
561 DSDV(14, 3,i)=(-4.5)+(36.)*z+(4.5)*x*x+(4.5)*y*y+(-0.405e2)*z*z+(-18.)*x*z+(-18.)*y*z;
562 DSDV(15, 3,i)=(-4.5)*x*x+(18.)*x*z;
563 DSDV(16, 3,i)=(-4.5)*y*y+(18.)*y*z;
564 }
565 #undef NUMSHAPES
566 #undef DIM
567 }
568
569
570 #undef V
571 #undef S
572 #undef DSDV
573
574

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26