/[escript]/branches/doubleplusgood/finley/src/ShapeFunctions.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/finley/src/ShapeFunctions.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4327 - (show annotations)
Wed Mar 20 05:09:11 2013 UTC (6 years, 7 months ago) by jfenwick
File size: 76033 byte(s)
some finley memory
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Finley: Shape functions */
20
21 /************************************************************************************/
22
23 #include "ShapeFunctions.h"
24 #include "esysUtils/mem.h"
25 #include "esysUtils/index.h"
26 #include <string.h>
27
28
29
30 Finley_ShapeFunctionInfo Finley_ShapeFunction_InfoList[]={
31 {Point1Shape, "Point1", 0, 1, 1, 1, Finley_Shape_Point1 } ,
32 {Line2Shape, "Line2", 1, 2, 1, 2, Finley_Shape_Line2 } ,
33 {Line3Shape, "Line3", 1, 3, 2, 2, Finley_Shape_Line3 },
34 {Line4Shape, "Line4", 1, 4, 3, 2, Finley_Shape_Line4 },
35 {Tri3Shape, "Tri3", 2, 3, 1, 3, Finley_Shape_Tri3 },
36 {Tri6Shape, "Tri6", 2, 6, 2, 3, Finley_Shape_Tri6 },
37 {Tri9Shape, "Tri9", 2, 9, 3, 3, Finley_Shape_Tri9 },
38 {Tri10Shape, "Tri10", 2, 10, 3, 3, Finley_Shape_Tri10, },
39 {Rec4Shape, "Rec4", 2, 4, 1, 4, Finley_Shape_Rec4, },
40 {Rec8Shape, "Rec8", 2, 8, 2, 4, Finley_Shape_Rec8, },
41 {Rec9Shape, "Rec9", 2, 9, 2, 4, Finley_Shape_Rec9, },
42 {Rec12Shape, "Rec12", 2, 12, 3, 4, Finley_Shape_Rec12, },
43 {Rec16Shape, "Rec16", 2, 16, 3, 4, Finley_Shape_Rec16, },
44 {Tet4Shape, "Tet4", 3, 4, 1, 4, Finley_Shape_Tet4, },
45 {Tet10Shape, "Tet10", 3, 10, 2, 4, Finley_Shape_Tet10, },
46 {Tet16Shape, "Tet16", 3, 16, 3, 4, Finley_Shape_Tet16, },
47 {Hex8Shape, "Hex8", 3, 8, 1, 8, Finley_Shape_Hex8, },
48 {Hex20Shape, "Hex20", 3, 20, 2, 8, Finley_Shape_Hex20, },
49 {Hex27Shape, "Hex27", 3, 27, 2, 8, Finley_Shape_Hex27, },
50 {Hex32Shape, "Hex32", 3, 32, 3, 8, Finley_Shape_Hex32, },
51 {NoShape, "NoType", 0, 1, 1, 1, Finley_Shape_Point1 }
52 };
53
54
55 /****************************************************************************************************************************************************
56
57 Creates an evaluation of the ShapeFunction on the given quadrature scheme.
58 if the spatial dimension of the scheme and the shape functions don't match
59
60 if QuadNodes==Null or QuadWeights==Null the shape functions method is used to generate a quadrature scheme with numQuadNodes
61 nodes. Otherwise it is assumed that a quadrature scheme is given on this array and a copy is created within the structure.
62
63 */
64 Finley_ShapeFunction* Finley_ShapeFunction_alloc(Finley_ShapeFunctionTypeId id,int numQuadDim, int numQuadNodes, double *QuadNodes, double *QuadWeights) {
65 Finley_ShapeFunction *out=NULL;
66 int numDim, numShapes, i, q;
67
68 numDim=Finley_ShapeFunction_InfoList[id].numDim;
69 numShapes=Finley_ShapeFunction_InfoList[id].numShapes;
70
71 if (numQuadDim>numDim) {
72 Finley_setError(VALUE_ERROR,"Finley_ShapeFunction_alloc: spatial dimension of quadrature scheme is bigger then spatial dimension of shape function.");
73 return NULL;
74 }
75
76 /* allocate the Finley_ShapeFunction to be returned: */
77
78 out=new Finley_ShapeFunction;
79 if (Finley_checkPtr(out)) return NULL;
80
81
82 out->Type=Finley_ShapeFunction_getInfo(id);
83 out->numQuadNodes=numQuadNodes;
84 out->QuadNodes=NULL;
85 out->QuadWeights=NULL;
86 out->S=NULL;
87 out->dSdv=NULL;
88 out->reference_counter=0;
89
90 /* allocate memory: */
91
92 out->QuadNodes=new double[numQuadNodes*numDim];
93 out->QuadWeights=new double[numQuadNodes];
94 out->S=new double[numShapes*numQuadNodes];
95 out->dSdv=new double[numShapes*numDim*numQuadNodes];
96 if ( Finley_checkPtr(out->QuadNodes) || Finley_checkPtr(out->QuadWeights) || Finley_checkPtr(out->S) || Finley_checkPtr(out->dSdv) ) {
97 Finley_ShapeFunction_dealloc(out);
98 return NULL;
99 }
100
101 /* set the quadrature nodes (missing values are filled with 0): */
102
103 for (q=0;q<numQuadNodes;q++) {
104 for (i=0;i<numQuadDim;i++) out->QuadNodes[INDEX2(i,q,numDim)]=QuadNodes[INDEX2(i,q,numQuadDim)];
105 for (i=numQuadDim;i<numDim;i++) out->QuadNodes[INDEX2(i,q,numDim)]=0;
106 out->QuadWeights[q]=QuadWeights[q];
107 }
108
109 /* eval shape functions on quadrature node: */
110
111 out->Type->getValues(numQuadNodes,out->QuadNodes,out->S,out->dSdv);
112
113 if (! Finley_noError()) {
114 Finley_ShapeFunction_dealloc(out);
115 return NULL;
116 }
117
118 /* all done: */
119 out->reference_counter=1;
120 return out;
121 }
122
123 Finley_ShapeFunction* Finley_ShapeFunction_reference(Finley_ShapeFunction* in) {
124 if (in!=NULL) ++(in->reference_counter);
125 return in;
126 }
127 /************************************************************************************/
128
129 void Finley_ShapeFunction_dealloc(Finley_ShapeFunction* in) {
130 if (in!=NULL) {
131 in->reference_counter--;
132 if (in->reference_counter<1) {
133 delete[] in->QuadNodes;
134 delete[] in->QuadWeights;
135 delete[] in->S;
136 delete[] in->dSdv;
137 delete in;
138 }
139 }
140 }
141
142 /************************************************************************************/
143
144 Finley_ShapeFunctionTypeId Finley_ShapeFunction_getTypeId(char* element_type)
145 {
146 int ptr=0;
147 Finley_ShapeFunctionTypeId out=NoShape;
148 while (Finley_ShapeFunction_InfoList[ptr].TypeId!=NoShape && out==NoShape) {
149 if (strcmp(element_type,Finley_ShapeFunction_InfoList[ptr].Name)==0) out=Finley_ShapeFunction_InfoList[ptr].TypeId;
150 ptr++;
151 }
152 return out;
153 }
154
155 Finley_ShapeFunctionInfo* Finley_ShapeFunction_getInfo(Finley_ShapeFunctionTypeId id)
156 {
157 int ptr=0;
158 Finley_ShapeFunctionInfo* out=NULL;
159 while (Finley_ShapeFunction_InfoList[ptr].TypeId!=NoShape && out==NULL) {
160 if (Finley_ShapeFunction_InfoList[ptr].TypeId==id) out=&(Finley_ShapeFunction_InfoList[ptr]);
161 ptr++;
162 }
163 if (out==NULL) {
164 Finley_setError(VALUE_ERROR,"Finley_ShapeFunctionInfo_getInfo: cannot find requested shape function");
165 }
166 return out;
167 }
168
169 /************************************************************************************/
170
171 #define V(_K_,_I_) v[INDEX2((_K_)-1,(_I_),DIM)]
172 #define S(_J_,_I_) s[S_INDEX((_J_)-1,(_I_),NUMSHAPES)]
173 #define DSDV(_J_,_K_,_I_) dsdv[DSDV_INDEX((_J_)-1,(_K_)-1,(_I_),NUMSHAPES,DIM)]
174
175 /************************************************************************************/
176
177 void Finley_Shape_Point1(int NumV,double* v,double* s,double* dsdv) {
178 #define NUMSHAPES 1
179 #define DIM 0
180 int i;
181 #pragma ivdep
182 for (i=0;i<NumV;i++) {
183 S(1,i)=1.;
184 }
185 #undef NUMSHAPES
186 #undef DIM
187 }
188
189 /************************************************************************************/
190
191 void Finley_Shape_Line2(int NumV,double* v,double* s,double* dsdv) {
192 #define NUMSHAPES 2
193 #define DIM 1
194 register double x;
195 int i;
196 #pragma ivdep
197 for (i=0;i<NumV;i++) {
198 x=V(1,i);
199 S(1,i)=1.-x;
200 S(2,i)= x;
201 DSDV(1,1,i)=-1.;
202 DSDV(2,1,i)= 1.;
203 }
204 #undef NUMSHAPES
205 #undef DIM
206 }
207
208 /************************************************************************************/
209
210 void Finley_Shape_Line3(int NumV,double* v,double* s,double* dsdv) {
211 #define NUMSHAPES 3
212 #define DIM 1
213 register double x;
214 int i;
215 #pragma ivdep
216 for (i=0;i<NumV;i++) {
217 x=V(1,i);
218 S(1,i)=(2.*x -1. )*(x -1.);
219 S(2,i)=(2.*x -1.)*x;
220 S(3,i)= 4.*x*(1. -x );
221 DSDV(1,1,i)= 4.*x -3.;
222 DSDV(2,1,i)= 4.*x -1.;
223 DSDV(3,1,i)=-8.*x+4.;
224 }
225 #undef NUMSHAPES
226 #undef DIM
227 }
228
229 /************************************************************************************/
230
231 void Finley_Shape_Line4(int NumV,double* v,double* s,double* dsdv) {
232 #define NUMSHAPES 4
233 #define DIM 1
234 register double x;
235 int i;
236 #pragma ivdep
237 for (i=0;i<NumV;i++) {
238 x=V(1,i);
239 S(1,i)=(10.)+(-5.5)*x+(9.)*x*x+(-4.5)*x*x*x ;
240 S(2,i)=(10.)*x+(-4.5)*x*x+(4.5)*x*x*x ;
241 S(3,i)=(9.)*x+(-22.5)*x*x+(13.5)*x*x*x ;
242 S(4,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x;
243 DSDV(1,1,i)=(-5.5)+(18.)*x+(-13.5)*x*x;
244 DSDV(2,1,i)=(10.)+(-9.)*x+(13.5)*x*x;
245 DSDV(3,1,i)=(9.)+(-45.)*x+(0.405e2)*x*x;
246 DSDV(4,1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x;
247 }
248 #undef NUMSHAPES
249 #undef DIM
250 }
251
252 /************************************************************************************/
253
254 void Finley_Shape_Tri3(int NumV,double* v,double* s,double* dsdv) {
255 #define NUMSHAPES 3
256 #define DIM 2
257 register double x,y;
258 int i;
259 #pragma ivdep
260 for (i=0;i<NumV;i++) {
261 x=V(1,i);
262 y=V(2,i);
263 S(1,i)=1.-x-y;
264 S(2,i)= x;
265 S(3,i)= y;
266 DSDV(1,1,i)=-1.;
267 DSDV(1,2,i)=-1.;
268 DSDV(2,1,i)= 1.;
269 DSDV(2,2,i)= 0.;
270 DSDV(3,1,i)= 0.;
271 DSDV(3,2,i)= 1.;
272 }
273 #undef NUMSHAPES
274 #undef DIM
275 }
276
277 /************************************************************************************/
278
279 void Finley_Shape_Tri6(int NumV,double* v,double* s,double* dsdv) {
280 #define NUMSHAPES 6
281 #define DIM 2
282 register double x,y;
283 int i;
284 #pragma ivdep
285 for (i=0;i<NumV;i++) {
286 x=V(1,i);
287 y=V(2,i);
288 S(1,i)= (1. -x -y)*(1. -2.*x -2.* y);
289 S(2,i)= x*(2.* x -1.);
290 S(3,i)= y*(2.* y -1.);
291 S(4,i)= (1. -x -y)*4.* x;
292 S(5,i)= 4.*x*y;
293 S(6,i)= (1. -x -y)*4.* y;
294 DSDV(1,1,i)= -3.+4.*x+4.*y;
295 DSDV(1,2,i)= -3.+4.*x+4.*y;
296 DSDV(2,1,i)= -1.+4.*x;
297 DSDV(2,2,i)= 0.;
298 DSDV(3,1,i)= 0.;
299 DSDV(3,2,i)= -1. +4.*y;
300 DSDV(4,1,i)= 4. -8.*x -4.*y;
301 DSDV(4,2,i)= -4.*x;
302 DSDV(5,1,i)= 4.*y;
303 DSDV(5,2,i)= 4.*x;
304 DSDV(6,1,i)= -4.*y;
305 DSDV(6,2,i)= 4. -4.*x -8.*y;
306 }
307 #undef NUMSHAPES
308 #undef DIM
309 }
310
311 /************************************************************************************/
312
313 void Finley_Shape_Tri9(int NumV,double* v,double* s,double* dsdv) {
314 #define NUMSHAPES 9
315 #define DIM 2
316 register double x,y;
317 int i;
318 #pragma ivdep
319 for (i=0;i<NumV;i++) {
320 x=V(1,i);
321 y=V(2,i);
322 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;
323 S(2,i)=(10.)*x+(-4.5)*x*x+(4.5)*x*x*x;
324 S(3,i)=(10.)*y+(-4.5)*y*y+(4.5)*y*y*y;
325 S(4,i)=(9.)*x+(-22.5)*x*x+(13.5)*x*x*x+(-9.)*x*y*y+(4.5)*x*x*y;
326 S(5,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x+(4.5)*x*y*y+(-9.)*x*x*y;
327 S(6,i)=(-4.5)*x*y*y+(9.)*x*x*y;
328 S(7,i)=(9.)*x*y*y+(-4.5)*x*x*y;
329 S(8,i)=(-4.5)*y+(18.)*y*y+(-13.5)*y*y*y+(-9.)*x*y*y+(4.5)*x*x*y;
330 S(9,i)=(9.)*y+(-22.5)*y*y+(13.5)*y*y*y+(4.5)*x*y*y+(-9.)*x*x*y;
331 DSDV(1, 1,i)=(-5.5)+(18.)*x+(-13.5)*x*x+(4.5)*y*y+(9.)*x*y;
332 DSDV(2, 1,i)=(10.)+(-9.)*x+(13.5)*x*x;
333 DSDV(3, 1,i)= 0.;
334 DSDV(4, 1,i)=(9.)+(-45.)*x+(0.405e2)*x*x+(-9.)*y*y+(9.)*x*y;
335 DSDV(5, 1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x+(4.5)*y*y+(-18.)*x*y;
336 DSDV(6, 1,i)=(-4.5)*y*y+(18.)*x*y;
337 DSDV(7, 1,i)=(9.)*y*y+(-9.)*x*y;
338 DSDV(8, 1,i)=(-9.)*y*y+(9.)*x*y;
339 DSDV(9, 1,i)=(4.5)*y*y+(-18.)*x*y;
340 DSDV(1, 2,i)=(-5.5)+(18.)*y+(-13.5)*y*y+(9.)*x*y+(4.5)*x*x;
341 DSDV(2, 2,i)= 0.;
342 DSDV(3, 2,i)=(10.)+(-9.)*y+(13.5)*y*y;
343 DSDV(4, 2,i)=(-18.)*x*y+(4.5)*x*x;
344 DSDV(5, 2,i)=(9.)*x*y+(-9.)*x*x;
345 DSDV(6, 2,i)=(-9.)*x*y+(9.)*x*x;
346 DSDV(7, 2,i)=(18.)*x*y+(-4.5)*x*x;
347 DSDV(8, 2,i)=(-4.5)+(36.)*y+(-0.405e2)*y*y+(-18.)*x*y+(4.5)*x*x;
348 DSDV(9, 2,i)=(9.)+(-45.)*y+(0.405e2)*y*y+(9.)*x*y+(-9.)*x*x;
349 }
350 #undef NUMSHAPES
351 #undef DIM
352 }
353
354 /************************************************************************************/
355
356 void Finley_Shape_Tri10(int NumV,double* v,double* s,double* dsdv) {
357 #define NUMSHAPES 10
358 #define DIM 2
359 register double x,y;
360 int i;
361 #pragma ivdep
362 for (i=0;i<NumV;i++) {
363 x=V(1,i);
364 y=V(2,i);
365 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;
366 S(2,i)=(10.)*x+(-4.5)*x*x+(4.5)*x*x*x;
367 S(3,i)=(10.)*y+(-4.5)*y*y+(4.5)*y*y*y;
368 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;
369 S(5,i)=(-4.5)*x+(18.)*x*x+(-13.5)*x*x*x+(-13.5)*x*x*y+(4.5)*x*y;
370 S(6,i)=(13.5)*x*x*y+(-4.5)*x*y;
371 S(7,i)=(13.5)*x*y*y+(-4.5)*x*y;
372 S(8,i)=(-4.5)*y+(18.)*y*y+(-13.5)*y*y*y+(-13.5)*x*y*y+(4.5)*x*y;
373 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;
374 S(10,i)=(-0.27e2)*x*y*y+(-0.27e2)*x*x*y+(0.27e2)*x*y;
375 DSDV(1, 1,i)=(-5.5)+(18.)*x+(-13.5)*x*x+(-13.5)*y*y+(-0.27e2)*x*y+(18.)*y;
376 DSDV(2, 1,i)=(10.)+(-9.)*x+(13.5)*x*x;
377 DSDV(3, 1,i)= 0.;
378 DSDV(4, 1,i)=(9.)+(-45.)*x+(0.405e2)*x*x+(13.5)*y*y+(0.54e2)*x*y+(-22.5)*y;
379 DSDV(5, 1,i)=(-4.5)+(36.)*x+(-0.405e2)*x*x+(-0.27e2)*x*y+(4.5)*y;
380 DSDV(6, 1,i)=(0.27e2)*x*y+(-4.5)*y;
381 DSDV(7, 1,i)=(13.5)*y*y+(-4.5)*y;
382 DSDV(8, 1,i)=(-13.5)*y*y+(4.5)*y;
383 DSDV(9, 1,i)=(0.27e2)*y*y+(0.27e2)*x*y+(-22.5)*y;
384 DSDV(10, 1,i)=(-0.27e2)*y*y+(-0.54e2)*x*y+(0.27e2)*y;
385 DSDV(1, 2,i)=(-5.5)+(18.)*y+(-13.5)*y*y+(-0.27e2)*x*y+(-13.5)*x*x+(18.)*x;
386 DSDV(2, 2,i)=0.;
387 DSDV(3, 2,i)=(10.)+(-9.)*y+(13.5)*y*y;
388 DSDV(4, 2,i)=(0.27e2)*x*y+(0.27e2)*x*x+(-22.5)*x;
389 DSDV(5, 2,i)=(-13.5)*x*x+(4.5)*x;
390 DSDV(6, 2,i)=(13.5)*x*x+(-4.5)*x;
391 DSDV(7, 2,i)=(0.27e2)*x*y+(-4.5)*x;
392 DSDV(8, 2,i)=(-4.5)+(36.)*y+(-0.405e2)*y*y+(-0.27e2)*x*y+(4.5)*x;
393 DSDV(9, 2,i)=(9.)+(-45.)*y+(0.405e2)*y*y+(0.54e2)*x*y+(13.5)*x*x+(-22.5)*x;
394 DSDV(10, 2,i)=(-0.54e2)*x*y+(-0.27e2)*x*x+(0.27e2)*x;
395 }
396 #undef NUMSHAPES
397 #undef DIM
398 }
399
400 /************************************************************************************/
401
402 void Finley_Shape_Rec4(int NumV,double* v,double* s,double* dsdv) {
403 #define NUMSHAPES 4
404 #define DIM 2
405 register double x,y;
406 int i;
407 #pragma ivdep
408 for (i=0;i<NumV;i++) {
409 x=V(1,i);
410 y=V(2,i);
411 S(1,i)=(1.-x)*(1.-y);
412 S(2,i)= x*(1.-y);
413 S(3,i)= x*y;
414 S(4,i)= (1.-x)*y;
415 DSDV(1,1,i)=y-1.;
416 DSDV(1,2,i)=x-1.;
417 DSDV(2,1,i)= 1.-y;
418 DSDV(2,2,i)=-x;
419 DSDV(3,1,i)= y;
420 DSDV(3,2,i)= x;
421 DSDV(4,1,i)=-y;
422 DSDV(4,2,i)= 1.-x;
423 }
424 #undef NUMSHAPES
425 #undef DIM
426 }
427
428 /************************************************************************************/
429
430 void Finley_Shape_Rec8(int NumV,double* v,double* s,double* dsdv) {
431 #define NUMSHAPES 8
432 #define DIM 2
433 register double x,y;
434 int i;
435 #pragma ivdep
436 for (i=0;i<NumV;i++) {
437 x=V(1,i);
438 y=V(2,i);
439 S(1,i)= 1. -3.*(x+y )+2.*x*x*(1. -y )+2.*y*y*(1. -x )+5.*x*y;
440 S(2,i)= x*(-1. -y+2.*x+2.*y*y -2.*x*y );
441 S(3,i)= x*y*(-3.+2.*(x+y ));
442 S(4,i)= y*(-1. -x+2.*y+2.*x*x -2.*x*y );
443 S(5,i)=4.*x*(1. -x -y+x* y);
444 S(6,i)= 4.*x*y*(1. -y );
445 S(7,i)= 4.*x*y*(1. -x );
446 S(8,i)=4.*y*(1. -x -y+x* y);
447 DSDV(1,1,i)=-3.+4.*x*(1. -y )+y*(5. -2.*y );
448 DSDV(1,2,i)=-3.+4.*y*(1. -x )+x*(5. -2.*x );
449 DSDV(2,1,i)=-1.+4.*x*(1. -y )+y*(-1.+2.*y );
450 DSDV(2,2,i)= x*(-1. -2.*x+4.*y);
451 DSDV(3,1,i)= y*(-3.+4.*x+2.*y);
452 DSDV(3,2,i)= x*(-3.+4.*y+2.*x);
453 DSDV(4,1,i)= y*(-1. -2.*y+4.*x);
454 DSDV(4,2,i)=-1.+4.*y*(1. -x )+x*(-1.+2.*x );
455 DSDV(5,1,i)= 4.*(1. -y )+8.*x*(y -1. );
456 DSDV(5,2,i)= 4.*x*(x -1.);
457 DSDV(6,1,i)= 4.*y*(1. -y );
458 DSDV(6,2,i)= 4.*x*(1. -2.*y );
459 DSDV(7,1,i)= 4.*y*(1. -2.*x );
460 DSDV(7,2,i)= 4.*x*(1. -x );
461 DSDV(8,1,i)= 4.*y*(y -1.);
462 DSDV(8,2,i)= 4.*(1. -x )+8.*y*(x -1. );
463 }
464 #undef NUMSHAPES
465 #undef DIM
466 }
467
468 /************************************************************************************/
469
470 void Finley_Shape_Rec9(int NumV,double* v,double* s,double* dsdv) {
471 #define NUMSHAPES 9
472 #define DIM 2
473 register double x,y;
474 int i;
475 #pragma ivdep
476 for (i=0;i<NumV;i++) {
477 x=V(1,i);
478 y=V(2,i);
479 S(1,i)= + 1.0 - 3.0*x + 2.0*x*x - 3.0*y + 9.0*x*y - 6.0*x*x*y + 2.0*y*y - 6.0*x*y*y + 4.0*x*x*y*y;
480 S(2,i)= - 1.0*x + 2.0*x*x + 3.0*x*y - 6.0*x*x*y - 2.0*x*y*y + 4.0*x*x*y*y;
481 S(3,i)= + 1.0*x*y - 2.0*x*x*y - 2.0*x*y*y + 4.0*x*x*y*y;
482 S(4,i)= - 1.0*y + 3.0*x*y - 2.0*x*x*y + 2.0*y*y - 6.0*x*y*y + 4.0*x*x*y*y;
483 S(5,i)= + 4.0*x - 4.0*x*x - 12.0*x*y + 12.0*x*x*y + 8.0*x*y*y - 8.0*x*x*y*y;
484 S(6,i)= - 4.0*x*y + 8.0*x*x*y + 4.0*x*y*y - 8.0*x*x*y*y;
485 S(7,i)= - 4.0*x*y + 4.0*x*x*y + 8.0*x*y*y - 8.0*x*x*y*y;
486 S(8,i)= + 4.0*y - 12.0*x*y + 8.0*x*x*y - 4.0*y*y + 12.0*x*y*y - 8.0*x*x*y*y;
487 S(9,i)= + 16.0*x*y - 16.0*x*x*y - 16.0*x*y*y + 16.0*x*x*y*y;
488 DSDV(1,1,i)= - 3.0 + 4.0*x + 9.0*y - 12.0*x*y - 6.0*y*y + 8.0*x*y*y;
489 DSDV(1,2,i)= - 3.0 + 9.0*x - 6.0*x*x + 4.0*y - 12.0*x*y + 8.0*x*x*y;
490 DSDV(2,1,i)= - 1.0 + 4.0*x + 3.0*y - 12.0*x*y - 2.0*y*y + 8.0*x*y*y;
491 DSDV(2,2,i)= + 3.0*x - 6.0*x*x - 4.0*x*y + 8.0*x*x*y;
492 DSDV(3,1,i)= + 1.0*y - 4.0*x*y - 2.0*y*y + 8.0*x*y*y;
493 DSDV(3,2,i)= + 1.0*x - 2.0*x*x - 4.0*x*y + 8.0*x*x*y;
494 DSDV(4,1,i)= + 3.0*y - 4.0*x*y - 6.0*y*y + 8.0*x*y*y;
495 DSDV(4,2,i)= - 1.0 + 3.0*x - 2.0*x*x + 4.0*y - 12.0*x*y + 8.0*x*x*y;
496 DSDV(5,1,i)= + 4.0 - 8.0*x - 12.0*y + 24.0*x*y + 8.0*y*y - 16.0*x*y*y;
497 DSDV(5,2,i)= - 12.0*x + 12.0*x*x + 16.0*x*y - 16.0*x*x*y;
498 DSDV(6,1,i)= - 4.0*y + 16.0*x*y + 4.0*y*y - 16.0*x*y*y;
499 DSDV(6,2,i)= - 4.0*x + 8.0*x*x + 8.0*x*y - 16.0*x*x*y;
500 DSDV(7,1,i)= - 4.0*y + 8.0*x*y + 8.0*y*y - 16.0*x*y*y;
501 DSDV(7,2,i)= - 4.0*x + 4.0*x*x + 16.0*x*y - 16.0*x*x*y;
502 DSDV(8,1,i)= - 12.0*y + 16.0*x*y + 12.0*y*y - 16.0*x*y*y;
503 DSDV(8,2,i)= + 4.0 - 12.0*x + 8.0*x*x - 8.0*y + 24.0*x*y - 16.0*x*x*y;
504 DSDV(9,1,i)= + 16.0*y - 32.0*x*y - 16.0*y*y + 32.0*x*y*y;
505 DSDV(9,2,i)= + 16.0*x - 16.0*x*x - 32.0*x*y + 32.0*x*x*y;
506 }
507 #undef NUMSHAPES
508 #undef DIM
509 }
510
511 /************************************************************************************/
512
513 void Finley_Shape_Rec12(int NumV,double* v,double* s,double* dsdv) {
514 #define NUMSHAPES 12
515 #define DIM 2
516 register double x,y;
517 int i;
518 #pragma ivdep
519 for (i=0;i<NumV;i++) {
520 x=V(1,i);
521 y=V(2,i);
522 S(1,i)=(10.)+(-5.5)*x+(10.)*x*y+(-5.5)*y+(9.)*x*x+(-4.5)*x*x*x+(-9.)*x*y*y+(4.5)*x*y*y*y+(4.5)*x*x*x*y+(-9.)*x*x*y+(-4.5)*y*y*y+(9.)*y*y;
523 S(2,i)=(10.)*x+(-5.5)*x*y+(-4.5)*x*x+(4.5)*x*x*x+(9.)*x*y*y+(-4.5)*x*y*y*y+(-4.5)*x*x*x*y+(4.5)*x*x*y;
524 S(3,i)=(10.)*x*y+(-4.5)*x*y*y+(4.5)*x*y*y*y+(4.5)*x*x*x*y+(-4.5)*x*x*y;
525 S(4,i)=(-5.5)*x*y+(10.)*y+(4.5)*x*y*y+(-4.5)*x*y*y*y+(-4.5)*x*x*x*y+(9.)*x*x*y+(4.5)*y*y*y+(-4.5)*y*y;
526 S(5,i)=(9.)*x+(-9.)*x*y+(-22.5)*x*x+(13.5)*x*x*x+(-13.5)*x*x*x*y+(22.5)*x*x*y;
527 S(6,i)=(-4.5)*x+(4.5)*x*y+(18.)*x*x+(-13.5)*x*x*x+(13.5)*x*x*x*y+(-18.)*x*x*y;
528 S(7,i)=(9.)*x*y+(-22.5)*x*y*y+(13.5)*x*y*y*y;
529 S(8,i)=(-4.5)*x*y+(18.)*x*y*y+(-13.5)*x*y*y*y;
530 S(9,i)=(-4.5)*x*y+(-13.5)*x*x*x*y+(18.)*x*x*y;
531 S(10,i)=(9.)*x*y+(13.5)*x*x*x*y+(-22.5)*x*x*y;
532 S(11,i)=(4.5)*x*y+(-4.5)*y+(-18.)*x*y*y+(13.5)*x*y*y*y+(-13.5)*y*y*y+(18.)*y*y;
533 S(12,i)=(-9.)*x*y+(9.)*y+(22.5)*x*y*y+(-13.5)*x*y*y*y+(13.5)*y*y*y+(-22.5)*y*y;
534 DSDV(1,1,i)=(-5.5)+(10.)*y+(18.)*x+(-13.5)*x*x+(-9.)*y*y+(4.5)*y*y*y+(13.5)*x*x*y+(-18.)*x*y;
535 DSDV(2,1,i)=(10.)+(-5.5)*y+(-9.)*x+(13.5)*x*x+(9.)*y*y+(-4.5)*y*y*y+(-13.5)*x*x*y+(9.)*x*y;
536 DSDV(3,1,i)=(10.)*y+(-4.5)*y*y+(4.5)*y*y*y+(13.5)*x*x*y+(-9.)*x*y;
537 DSDV(4,1,i)=(-5.5)*y+(4.5)*y*y+(-4.5)*y*y*y+(-13.5)*x*x*y+(18.)*x*y;
538 DSDV(5,1,i)=(9.)+(-9.)*y+(-45.)*x+(0.405e2)*x*x+(-0.405e2)*x*x*y+(45.)*x*y;
539 DSDV(6,1,i)=(-4.5)+(4.5)*y+(36.)*x+(-0.405e2)*x*x+(0.405e2)*x*x*y+(-36.)*x*y;
540 DSDV(7,1,i)=(9.)*y+(-22.5)*y*y+(13.5)*y*y*y;
541 DSDV(8,1,i)=(-4.5)*y+(18.)*y*y+(-13.5)*y*y*y;
542 DSDV(9,1,i)=(-4.5)*y+(-0.405e2)*x*x*y+(36.)*x*y;
543 DSDV(10,1,i)=(9.)*y+(0.405e2)*x*x*y+(-45.)*x*y;
544 DSDV(11,1,i)=(4.5)*y+(-18.)*y*y+(13.5)*y*y*y;
545 DSDV(12,1,i)=(-9.)*y+(22.5)*y*y+(-13.5)*y*y*y;
546 DSDV(1,2,i)=(10.)*x+(-5.5)+(-18.)*x*y+(13.5)*x*y*y+(4.5)*x*x*x+(-9.)*x*x+(-13.5)*y*y+(18.)*y;
547 DSDV(2,2,i)=(-5.5)*x+(18.)*x*y+(-13.5)*x*y*y+(-4.5)*x*x*x+(4.5)*x*x;
548 DSDV(3,2,i)=(10.)*x+(-9.)*x*y+(13.5)*x*y*y+(4.5)*x*x*x+(-4.5)*x*x;
549 DSDV(4,2,i)=(-5.5)*x+(10.)+(9.)*x*y+(-13.5)*x*y*y+(-4.5)*x*x*x+(9.)*x*x+(13.5)*y*y+(-9.)*y;
550 DSDV(5,2,i)=(-9.)*x+(-13.5)*x*x*x+(22.5)*x*x;
551 DSDV(6,2,i)=(4.5)*x+(13.5)*x*x*x+(-18.)*x*x;
552 DSDV(7,2,i)=(9.)*x+(-45.)*x*y+(0.405e2)*x*y*y;
553 DSDV(8,2,i)=(-4.5)*x+(36.)*x*y+(-0.405e2)*x*y*y;
554 DSDV(9,2,i)=(-4.5)*x+(-13.5)*x*x*x+(18.)*x*x;
555 DSDV(10,2,i)=(9.)*x+(13.5)*x*x*x+(-22.5)*x*x;
556 DSDV(11,2,i)=(4.5)*x+(-4.5)+(-36.)*x*y+(0.405e2)*x*y*y+(-0.405e2)*y*y+(36.)*y;
557 DSDV(12,2,i)=(-9.)*x+(9.)+(45.)*x*y+(-0.405e2)*x*y*y+(0.405e2)*y*y+(-45.)*y;
558 }
559 #undef NUMSHAPES
560 #undef DIM
561 }
562
563 /************************************************************************************/
564
565 void Finley_Shape_Rec16(int NumV,double* v,double* s,double* dsdv) {
566 #define NUMSHAPES 16
567 #define DIM 2
568 register double x,y;
569 int i;
570 #pragma ivdep
571 for (i=0;i<NumV;i++) {
572 x=V(1,i);
573 y=V(2,i);
574 S(1,i)=(10.)+(-5.5)*x+(0.3025e2)*x*y+(-5.5)*y+(9.)*x*x+(-4.5)*x*x*x+(-0.495e2)*x*y*y+(0.2475e2)*x*y*y*y+(0.2475e2)*x*x*x*y+(-0.495e2)*x*x*y+(-4.5)*y*y*y+(9.)*y*y+(0.81e2)*x*x*y*y+(-0.405e2)*x*x*x*y*y+(0.2025e2)*x*x*x*y*y*y+(-0.405e2)*x*x*y*y*y;
575 S(2,i)=(10.)*x+(-5.5)*x*y+(-4.5)*x*x+(4.5)*x*x*x+(9.)*x*y*y+(-4.5)*x*y*y*y+(-0.2475e2)*x*x*x*y+(0.2475e2)*x*x*y+(-0.405e2)*x*x*y*y+(0.405e2)*x*x*x*y*y+(-0.2025e2)*x*x*x*y*y*y+(0.2025e2)*x*x*y*y*y;
576 S(3,i)=(10.)*x*y+(-4.5)*x*y*y+(4.5)*x*y*y*y+(4.5)*x*x*x*y+(-4.5)*x*x*y+(0.2025e2)*x*x*y*y+(-0.2025e2)*x*x*x*y*y+(0.2025e2)*x*x*x*y*y*y+(-0.2025e2)*x*x*y*y*y;
577 S(4,i)=(-5.5)*x*y+(10.)*y+(0.2475e2)*x*y*y+(-0.2475e2)*x*y*y*y+(-4.5)*x*x*x*y+(9.)*x*x*y+(4.5)*y*y*y+(-4.5)*y*y+(-0.405e2)*x*x*y*y+(0.2025e2)*x*x*x*y*y+(-0.2025e2)*x*x*x*y*y*y+(0.405e2)*x*x*y*y*y;
578 S(5,i)=(9.)*x+(-0.495e2)*x*y+(-22.5)*x*x+(13.5)*x*x*x+(0.81e2)*x*y*y+(-0.405e2)*x*y*y*y+(-0.7425e2)*x*x*x*y+(0.12375e3)*x*x*y+(-0.2025e3)*x*x*y*y+(0.1215e3)*x*x*x*y*y+(-0.6075e2)*x*x*x*y*y*y+(0.10125e3)*x*x*y*y*y;
579 S(6,i)=(-4.5)*x+(0.2475e2)*x*y+(18.)*x*x+(-13.5)*x*x*x+(-0.405e2)*x*y*y+(0.2025e2)*x*y*y*y+(0.7425e2)*x*x*x*y+(-0.99e2)*x*x*y+(0.162e3)*x*x*y*y+(-0.1215e3)*x*x*x*y*y+(0.6075e2)*x*x*x*y*y*y+(-0.81e2)*x*x*y*y*y;
580 S(7,i)=(9.)*x*y+(-22.5)*x*y*y+(13.5)*x*y*y*y+(0.405e2)*x*x*x*y+(-0.405e2)*x*x*y+(0.10125e3)*x*x*y*y+(-0.10125e3)*x*x*x*y*y+(0.6075e2)*x*x*x*y*y*y+(-0.6075e2)*x*x*y*y*y;
581 S(8,i)=(-4.5)*x*y+(18.)*x*y*y+(-13.5)*x*y*y*y+(-0.2025e2)*x*x*x*y+(0.2025e2)*x*x*y+(-0.81e2)*x*x*y*y+(0.81e2)*x*x*x*y*y+(-0.6075e2)*x*x*x*y*y*y+(0.6075e2)*x*x*y*y*y;
582 S(9,i)=(-4.5)*x*y+(0.2025e2)*x*y*y+(-0.2025e2)*x*y*y*y+(-13.5)*x*x*x*y+(18.)*x*x*y+(-0.81e2)*x*x*y*y+(0.6075e2)*x*x*x*y*y+(-0.6075e2)*x*x*x*y*y*y+(0.81e2)*x*x*y*y*y;
583 S(10,i)=(9.)*x*y+(-0.405e2)*x*y*y+(0.405e2)*x*y*y*y+(13.5)*x*x*x*y+(-22.5)*x*x*y+(0.10125e3)*x*x*y*y+(-0.6075e2)*x*x*x*y*y+(0.6075e2)*x*x*x*y*y*y+(-0.10125e3)*x*x*y*y*y;
584 S(11,i)=(0.2475e2)*x*y+(-4.5)*y+(-0.99e2)*x*y*y+(0.7425e2)*x*y*y*y+(0.2025e2)*x*x*x*y+(-0.405e2)*x*x*y+(-13.5)*y*y*y+(18.)*y*y+(0.162e3)*x*x*y*y+(-0.81e2)*x*x*x*y*y+(0.6075e2)*x*x*x*y*y*y+(-0.1215e3)*x*x*y*y*y;
585 S(12,i)=(-0.495e2)*x*y+(9.)*y+(0.12375e3)*x*y*y+(-0.7425e2)*x*y*y*y+(-0.405e2)*x*x*x*y+(0.81e2)*x*x*y+(13.5)*y*y*y+(-22.5)*y*y+(-0.2025e3)*x*x*y*y+(0.10125e3)*x*x*x*y*y+(-0.6075e2)*x*x*x*y*y*y+(0.1215e3)*x*x*y*y*y;
586 S(13,i)=(0.81e2)*x*y+(-0.2025e3)*x*y*y+(0.1215e3)*x*y*y*y+(0.1215e3)*x*x*x*y+(-0.2025e3)*x*x*y+(0.50625e3)*x*x*y*y+(-0.30375e3)*x*x*x*y*y+(0.18225e3)*x*x*x*y*y*y+(-0.30375e3)*x*x*y*y*y;
587 S(14,i)=(-0.405e2)*x*y+(0.10125e3)*x*y*y+(-0.6075e2)*x*y*y*y+(-0.1215e3)*x*x*x*y+(0.162e3)*x*x*y+(-0.405e3)*x*x*y*y+(0.30375e3)*x*x*x*y*y+(-0.18225e3)*x*x*x*y*y*y+(0.243e3)*x*x*y*y*y;
588 S(15,i)=(0.2025e2)*x*y+(-0.81e2)*x*y*y+(0.6075e2)*x*y*y*y+(0.6075e2)*x*x*x*y+(-0.81e2)*x*x*y+(0.324e3)*x*x*y*y+(-0.243e3)*x*x*x*y*y+(0.18225e3)*x*x*x*y*y*y+(-0.243e3)*x*x*y*y*y;
589 S(16,i)=(-0.405e2)*x*y+(0.162e3)*x*y*y+(-0.1215e3)*x*y*y*y+(-0.6075e2)*x*x*x*y+(0.10125e3)*x*x*y+(-0.405e3)*x*x*y*y+(0.243e3)*x*x*x*y*y+(-0.18225e3)*x*x*x*y*y*y+(0.30375e3)*x*x*y*y*y;
590 DSDV(1, 1,i)=(-5.5)+(0.3025e2)*y+(18.)*x+(-13.5)*x*x+(-0.495e2)*y*y+(0.2475e2)*y*y*y+(0.7425e2)*x*x*y+(-0.99e2)*x*y+(0.162e3)*x*y*y+(-0.1215e3)*x*x*y*y+(0.6075e2)*x*x*y*y*y+(-0.81e2)*x*y*y*y;
591 DSDV(2, 1,i)=(10.)+(-5.5)*y+(-9.)*x+(13.5)*x*x+(9.)*y*y+(-4.5)*y*y*y+(-0.7425e2)*x*x*y+(0.495e2)*x*y+(-0.81e2)*x*y*y+(0.1215e3)*x*x*y*y+(-0.6075e2)*x*x*y*y*y+(0.405e2)*x*y*y*y;
592 DSDV(3, 1,i)=(10.)*y+(-4.5)*y*y+(4.5)*y*y*y+(13.5)*x*x*y+(-9.)*x*y+(0.405e2)*x*y*y+(-0.6075e2)*x*x*y*y+(0.6075e2)*x*x*y*y*y+(-0.405e2)*x*y*y*y;
593 DSDV(4, 1,i)=(-5.5)*y+(0.2475e2)*y*y+(-0.2475e2)*y*y*y+(-13.5)*x*x*y+(18.)*x*y+(-0.81e2)*x*y*y+(0.6075e2)*x*x*y*y+(-0.6075e2)*x*x*y*y*y+(0.81e2)*x*y*y*y;
594 DSDV(5, 1,i)=(9.)+(-0.495e2)*y+(-45.)*x+(0.405e2)*x*x+(0.81e2)*y*y+(-0.405e2)*y*y*y+(-0.22275e3)*x*x*y+(0.2475e3)*x*y+(-0.405e3)*x*y*y+(0.3645e3)*x*x*y*y+(-0.18225e3)*x*x*y*y*y+(0.2025e3)*x*y*y*y;
595 DSDV(6, 1,i)=(-4.5)+(0.2475e2)*y+(36.)*x+(-0.405e2)*x*x+(-0.405e2)*y*y+(0.2025e2)*y*y*y+(0.22275e3)*x*x*y+(-0.198e3)*x*y+(0.324e3)*x*y*y+(-0.3645e3)*x*x*y*y+(0.18225e3)*x*x*y*y*y+(-0.162e3)*x*y*y*y;
596 DSDV(7, 1,i)=(9.)*y+(-22.5)*y*y+(13.5)*y*y*y+(0.1215e3)*x*x*y+(-0.81e2)*x*y+(0.2025e3)*x*y*y+(-0.30375e3)*x*x*y*y+(0.18225e3)*x*x*y*y*y+(-0.1215e3)*x*y*y*y;
597 DSDV(8, 1,i)=(-4.5)*y+(18.)*y*y+(-13.5)*y*y*y+(-0.6075e2)*x*x*y+(0.405e2)*x*y+(-0.162e3)*x*y*y+(0.243e3)*x*x*y*y+(-0.18225e3)*x*x*y*y*y+(0.1215e3)*x*y*y*y;
598 DSDV(9, 1,i)=(-4.5)*y+(0.2025e2)*y*y+(-0.2025e2)*y*y*y+(-0.405e2)*x*x*y+(36.)*x*y+(-0.162e3)*x*y*y+(0.18225e3)*x*x*y*y+(-0.18225e3)*x*x*y*y*y+(0.162e3)*x*y*y*y;
599 DSDV(10, 1,i)=(9.)*y+(-0.405e2)*y*y+(0.405e2)*y*y*y+(0.405e2)*x*x*y+(-45.)*x*y+(0.2025e3)*x*y*y+(-0.18225e3)*x*x*y*y+(0.18225e3)*x*x*y*y*y+(-0.2025e3)*x*y*y*y;
600 DSDV(11, 1,i)=(0.2475e2)*y+(-0.99e2)*y*y+(0.7425e2)*y*y*y+(0.6075e2)*x*x*y+(-0.81e2)*x*y+(0.324e3)*x*y*y+(-0.243e3)*x*x*y*y+(0.18225e3)*x*x*y*y*y+(-0.243e3)*x*y*y*y;
601 DSDV(12, 1,i)=(-0.495e2)*y+(0.12375e3)*y*y+(-0.7425e2)*y*y*y+(-0.1215e3)*x*x*y+(0.162e3)*x*y+(-0.405e3)*x*y*y+(0.30375e3)*x*x*y*y+(-0.18225e3)*x*x*y*y*y+(0.243e3)*x*y*y*y;
602 DSDV(13, 1,i)=(0.81e2)*y+(-0.2025e3)*y*y+(0.1215e3)*y*y*y+(0.3645e3)*x*x*y+(-0.405e3)*x*y+(0.10125e4)*x*y*y+(-0.91125e3)*x*x*y*y+(0.54675e3)*x*x*y*y*y+(-0.6075e3)*x*y*y*y;
603 DSDV(14, 1,i)=(-0.405e2)*y+(0.10125e3)*y*y+(-0.6075e2)*y*y*y+(-0.3645e3)*x*x*y+(0.324e3)*x*y+(-0.81e3)*x*y*y+(0.91125e3)*x*x*y*y+(-0.54675e3)*x*x*y*y*y+(0.486e3)*x*y*y*y;
604 DSDV(15, 1,i)=(0.2025e2)*y+(-0.81e2)*y*y+(0.6075e2)*y*y*y+(0.18225e3)*x*x*y+(-0.162e3)*x*y+(0.648e3)*x*y*y+(-0.729e3)*x*x*y*y+(0.54675e3)*x*x*y*y*y+(-0.486e3)*x*y*y*y;
605 DSDV(16, 1,i)=(-0.405e2)*y+(0.162e3)*y*y+(-0.1215e3)*y*y*y+(-0.18225e3)*x*x*y+(0.2025e3)*x*y+(-0.81e3)*x*y*y+(0.729e3)*x*x*y*y+(-0.54675e3)*x*x*y*y*y+(0.6075e3)*x*y*y*y;
606 DSDV(1, 2,i)=(0.3025e2)*x+(-5.5)+(-0.99e2)*x*y+(0.7425e2)*x*y*y+(0.2475e2)*x*x*x+(-0.495e2)*x*x+(-13.5)*y*y+(18.)*y+(0.162e3)*x*x*y+(-0.81e2)*x*x*x*y+(0.6075e2)*x*x*x*y*y+(-0.1215e3)*x*x*y*y;
607 DSDV(2, 2,i)=(-5.5)*x+(18.)*x*y+(-13.5)*x*y*y+(-0.2475e2)*x*x*x+(0.2475e2)*x*x+(-0.81e2)*x*x*y+(0.81e2)*x*x*x*y+(-0.6075e2)*x*x*x*y*y+(0.6075e2)*x*x*y*y;
608 DSDV(3, 2,i)=(10.)*x+(-9.)*x*y+(13.5)*x*y*y+(4.5)*x*x*x+(-4.5)*x*x+(0.405e2)*x*x*y+(-0.405e2)*x*x*x*y+(0.6075e2)*x*x*x*y*y+(-0.6075e2)*x*x*y*y;
609 DSDV(4, 2,i)=(-5.5)*x+(10.)+(0.495e2)*x*y+(-0.7425e2)*x*y*y+(-4.5)*x*x*x+(9.)*x*x+(13.5)*y*y+(-9.)*y+(-0.81e2)*x*x*y+(0.405e2)*x*x*x*y+(-0.6075e2)*x*x*x*y*y+(0.1215e3)*x*x*y*y;
610 DSDV(5, 2,i)=(-0.495e2)*x+(0.162e3)*x*y+(-0.1215e3)*x*y*y+(-0.7425e2)*x*x*x+(0.12375e3)*x*x+(-0.405e3)*x*x*y+(0.243e3)*x*x*x*y+(-0.18225e3)*x*x*x*y*y+(0.30375e3)*x*x*y*y;
611 DSDV(6, 2,i)=(0.2475e2)*x+(-0.81e2)*x*y+(0.6075e2)*x*y*y+(0.7425e2)*x*x*x+(-0.99e2)*x*x+(0.324e3)*x*x*y+(-0.243e3)*x*x*x*y+(0.18225e3)*x*x*x*y*y+(-0.243e3)*x*x*y*y;
612 DSDV(7, 2,i)=(9.)*x+(-45.)*x*y+(0.405e2)*x*y*y+(0.405e2)*x*x*x+(-.405e2)*x*x+(0.2025e3)*x*x*y+(-0.2025e3)*x*x*x*y+(0.18225e3)*x*x*x*y*y+(-0.18225e3)*x*x*y*y;
613 DSDV(8, 2,i)=(-4.5)*x+(36.)*x*y+(-0.405e2)*x*y*y+(-0.2025e2)*x*x*x+(0.2025e2)*x*x+(-0.162e3)*x*x*y+(0.162e3)*x*x*x*y+(-0.18225e3)*x*x*x*y*y+(0.18225e3)*x*x*y*y;
614 DSDV(9, 2,i)=(-4.5)*x+(0.405e2)*x*y+(-0.6075e2)*x*y*y+(-13.5)*x*x*x+(18.)*x*x+(-0.162e3)*x*x*y+(0.1215e3)*x*x*x*y+(-0.18225e3)*x*x*x*y*y+(0.243e3)*x*x*y*y;
615 DSDV(10, 2,i)=(9.)*x+(-0.81e2)*x*y+(0.1215e3)*x*y*y+(13.5)*x*x*x+(-22.5)*x*x+(0.2025e3)*x*x*y+(-0.1215e3)*x*x*x*y+(0.18225e3)*x*x*x*y*y+(-0.30375e3)*x*x*y*y;
616 DSDV(11, 2,i)=(0.2475e2)*x+(-4.5)+(-0.198e3)*x*y+(0.22275e3)*x*y*y+(0.2025e2)*x*x*x+(-0.405e2)*x*x+(-0.405e2)*y*y+(36.)*y+(0.324e3)*x*x*y+(-0.162e3)*x*x*x*y+(0.18225e3)*x*x*x*y*y+(-0.3645e3)*x*x*y*y;
617 DSDV(12, 2,i)=(-0.495e2)*x+(9.)+(0.2475e3)*x*y+(-0.22275e3)*x*y*y+(-0.405e2)*x*x*x+(0.81e2)*x*x+(0.405e2)*y*y+(-45.)*y+(-0.405e3)*x*x*y+(0.2025e3)*x*x*x*y+(-0.18225e3)*x*x*x*y*y+(0.3645e3)*x*x*y*y;
618 DSDV(13, 2,i)=(0.81e2)*x+(-0.405e3)*x*y+(0.3645e3)*x*y*y+(0.1215e3)*x*x*x+(-0.2025e3)*x*x+(0.10125e4)*x*x*y+(-0.6075e3)*x*x*x*y+(0.54675e3)*x*x*x*y*y+(-0.91125e3)*x*x*y*y;
619 DSDV(14, 2,i)=(-0.405e2)*x+(0.2025e3)*x*y+(-0.18225e3)*x*y*y+(-0.1215e3)*x*x*x+(0.162e3)*x*x+(-0.81e3)*x*x*y+(0.6075e3)*x*x*x*y+(-0.54675e3)*x*x*x*y*y+(0.729e3)*x*x*y*y;
620 DSDV(15, 2,i)=(0.2025e2)*x+(-0.162e3)*x*y+(0.18225e3)*x*y*y+(0.6075e2)*x*x*x+(-0.81e2)*x*x+(0.648e3)*x*x*y+(-0.486e3)*x*x*x*y+(0.54675e3)*x*x*x*y*y+(-0.729e3)*x*x*y*y;
621 DSDV(16, 2,i)=(-0.405e2)*x+(0.324e3)*x*y+(-0.3645e3)*x*y*y+(-0.6075e2)*x*x*x+(0.10125e3)*x*x+(-0.81e3)*x*x*y+(0.486e3)*x*x*x*y+(-0.54675e3)*x*x*x*y*y+(0.91125e3)*x*x*y*y;
622 }
623 #undef NUMSHAPES
624 #undef DIM
625 }
626
627 /************************************************************************************/
628
629 void Finley_Shape_Tet4(int NumV,double* v,double* s,double* dsdv) {
630 #define NUMSHAPES 4
631 #define DIM 3
632 register double x,y,z;
633 int i;
634 #pragma ivdep
635 for (i=0;i<NumV;i++) {
636 x=V(1,i);
637 y=V(2,i);
638 z=V(3,i);
639 S(1,i)=1.-x-y-z;
640 S(2,i)=x;
641 S(3,i)=y;
642 S(4,i)=z;
643 DSDV(1,1,i)=-1.;
644 DSDV(1,2,i)=-1.;
645 DSDV(1,3,i)=-1.;
646 DSDV(2,1,i)= 1.;
647 DSDV(2,2,i)= 0.;
648 DSDV(2,3,i)= 0.;
649 DSDV(3,1,i)= 0.;
650 DSDV(3,2,i)= 1.;
651 DSDV(3,3,i)= 0.;
652 DSDV(4,1,i)= 0.;
653 DSDV(4,2,i)= 0.;
654 DSDV(4,3,i)= 1.;
655 }
656 #undef NUMSHAPES
657 #undef DIM
658 }
659
660 /************************************************************************************/
661
662 void Finley_Shape_Tet10(int NumV,double* v,double* s,double* dsdv) {
663 #define NUMSHAPES 10
664 #define DIM 3
665 register double x,y,z;
666 int i;
667 #pragma ivdep
668 for (i=0;i<NumV;i++) {
669 x=V(1,i);
670 y=V(2,i);
671 z=V(3,i);
672 S(1,i) = (1.-x-y-z)*(1.-2.*x-2.*y-2.*z);
673 S(2,i) = x*(2.*x-1.);
674 S(3,i) = y*(2.*y-1.);
675 S(4,i) = z*(2.*z-1.);
676 S(5,i) = (1.-x-y-z)*4.*x;
677 S(6,i) = 4.*x*y;
678 S(7,i) = (1.-x-y-z)*4.*y;
679 S(8,i) = (1.-x-y-z)*4.*z;
680 S(9,i) = 4.*x*z;
681 S(10,i)= 4.*y*z;
682
683 DSDV(1,1,i)= -3.+4.*x+4.*y+4.*z;
684 DSDV(1,2,i)= -3.+4.*x+4.*y+4.*z;
685 DSDV(1,3,i)= -3.+4.*x+4.*y+4.*z;
686
687
688 DSDV(2,1,i)= -1.+4.*x;
689 DSDV(2,2,i)= 0.;
690 DSDV(2,3,i)= 0.;
691
692 DSDV(3,1,i)= 0.;
693 DSDV(3,2,i)= -1. +4.*y;
694 DSDV(3,3,i)= 0.;
695
696 DSDV(4,1,i)= 0.;
697 DSDV(4,2,i)= 0.;
698 DSDV(4,3,i)= -1. +4.*z;
699
700 DSDV(5,1,i)= 4. -8.*x -4.*y -4.*z;
701 DSDV(5,2,i)= -4.*x;
702 DSDV(5,3,i)= -4.*x;
703
704 DSDV(6,1,i)= 4.*y;
705 DSDV(6,2,i)= 4.*x;
706 DSDV(6,3,i)= 0.;
707
708 DSDV(7,1,i)= -4.*y;
709 DSDV(7,2,i)= 4. -4.*x -8.*y -4.*z;
710 DSDV(7,3,i)= -4.*y;
711
712 DSDV(8,1,i)= -4.*z;
713 DSDV(8,2,i)= -4.*z;
714 DSDV(8,3,i)= 4. -4.*x -4.*y -8.*z;
715
716 DSDV(9,1,i)= 4.*z;
717 DSDV(9,2,i)= 0.;
718 DSDV(9,3,i)= 4.*x;
719
720 DSDV(10,1,i)= 0.;
721 DSDV(10,2,i)= 4.*z;
722 DSDV(10,3,i)= 4.*y;
723 }
724 #undef NUMSHAPES
725 #undef DIM
726 }
727
728 /************************************************************************************/
729
730 void Finley_Shape_Tet16(int NumV,double* v,double* s,double* dsdv) {
731 #define NUMSHAPES 16
732 #define DIM 3
733 register double x,y,z;
734 int i;
735 #pragma ivdep
736 for (i=0;i<NumV;i++) {
737 x=V(1,i);
738 y=V(2,i);
739 z=V(3,i);
740 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;
741 S(2,i)=(1.e0)*x+(-4.5)*x*x+(4.5)*x*x*x;
742 S(3,i)=(1.e0)*y+(4.5)*y*y*y+(-4.5)*y*y;
743 S(4,i)=(1.e0)*z+(-4.5)*z*z+(4.5)*z*z*z;
744 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;
745 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;
746 S(7,i)=(9.)*x*x*y+(-4.5)*x*y*y;
747 S(8,i)=(-4.5)*x*x*y+(9.)*x*y*y;
748 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;
749 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;
750 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;
751 S(12,i)=(9.)*x*x*z+(-4.5)*x*z*z;
752 S(13,i)=(9.)*y*y*z+(-4.5)*y*z*z;
753 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;
754 S(15,i)=(-4.5)*x*x*z+(9.)*x*z*z;
755 S(16,i)=(-4.5)*y*y*z+(9.)*y*z*z;
756 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;
757 DSDV(2, 1,i)=(1.e0)+(-9.)*x+(13.5)*x*x;
758 DSDV(3, 1,i)= 0.;
759 DSDV(4, 1,i)= 0.;
760 DSDV(5, 1,i)=(9.)+(-45.)*x+(0.405e2)*x*x+(9.)*x*y+(-9.)*y*y+(9.)*x*z+(-9.)*z*z;
761 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;
762 DSDV(7, 1,i)=(18.)*x*y+(-4.5)*y*y;
763 DSDV(8, 1,i)=(-9.)*x*y+(9.)*y*y;
764 DSDV(9, 1,i)=(9.)*x*y+(-9.)*y*y;
765 DSDV(10, 1,i)=(-18.)*x*y+(4.5)*y*y;
766 DSDV(11, 1,i)=(-18.)*x*z+(4.5)*z*z;
767 DSDV(12, 1,i)=(18.)*x*z+(-4.5)*z*z;
768 DSDV(13, 1,i)=0.;
769 DSDV(14, 1,i)=(9.)*x*z+(-9.)*z*z;
770 DSDV(15, 1,i)=(-9.)*x*z+(9.)*z*z;
771 DSDV(16, 1,i)=0.;
772 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;
773 DSDV(2, 2,i)=0.;
774 DSDV(3, 2,i)=(1.e0)+(13.5)*y*y+(-9.)*y;
775 DSDV(4, 2,i)=0.;
776 DSDV(5, 2,i)=(4.5)*x*x+(-18.)*x*y;
777 DSDV(6, 2,i)=(-9.)*x*x+(9.)*x*y;
778 DSDV(7, 2,i)=(9.)*x*x+(-9.)*x*y;
779 DSDV(8, 2,i)=(-4.5)*x*x+(18.)*x*y;
780 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;
781 DSDV(10, 2,i)=(9.)+(-9.)*x*x+(9.)*x*y+(0.405e2)*y*y+(-45.)*y+(9.)*y*z+(-9.)*z*z;
782 DSDV(11, 2,i)=(-18.)*y*z+(4.5)*z*z;
783 DSDV(12, 2,i)=0.;
784 DSDV(13, 2,i)=(18.)*y*z+(-4.5)*z*z;
785 DSDV(14, 2,i)=(9.)*y*z+(-9.)*z*z;
786 DSDV(15, 2,i)=0.;
787 DSDV(16, 2,i)=(-9.)*y*z+(9.)*z*z;
788 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;
789 DSDV(2, 3,i)= 0.;
790 DSDV(3, 3,i)= 0.;
791 DSDV(4, 3,i)=(1.e0)+(-9.)*z+(13.5)*z*z;
792 DSDV(5, 3,i)=(4.5)*x*x+(-18.)*x*z;
793 DSDV(6, 3,i)=(-9.)*x*x+(9.)*x*z;
794 DSDV(7, 3,i)= 0.;
795 DSDV(8, 3,i)= 0.;
796 DSDV(9, 3,i)=(-9.)*y*y+(9.)*y*z;
797 DSDV(10, 3,i)=(4.5)*y*y+(-18.)*y*z;
798 DSDV(11, 3,i)=(9.)+(-45.)*z+(-9.)*x*x+(-9.)*y*y+(0.405e2)*z*z+(.9e1)*x*z+(9.)*y*z;
799 DSDV(12, 3,i)=(9.)*x*x+(-9.)*x*z;
800 DSDV(13, 3,i)=(9.)*y*y+(-9.)*y*z;
801 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;
802 DSDV(15, 3,i)=(-4.5)*x*x+(18.)*x*z;
803 DSDV(16, 3,i)=(-4.5)*y*y+(18.)*y*z;
804 }
805 #undef NUMSHAPES
806 #undef DIM
807 }
808
809 /************************************************************************************/
810
811 void Finley_Shape_Hex8(int NumV,double* v,double* s,double* dsdv) {
812 #define NUMSHAPES 8
813 #define DIM 3
814 register double x,y,z;
815 int i;
816 #pragma ivdep
817 for (i=0;i<NumV;i++) {
818 x=V(1,i);
819 y=V(2,i);
820 z=V(3,i);
821 S(1,i)=(1.-x)*(1.-y)*(1.-z);
822 S(2,i)= x*(1.-z)*(1.-y);
823 S(3,i)= x*(1.-z)*y;
824 S(4,i)= (1.-z)*(1.-x)*y;
825 S(5,i)= (1.-x)*z*(1.-y);
826 S(6,i)= x*z*(1.-y);
827 S(7,i)= x*y*z;
828 S(8,i)= y*z*(1.-x);
829 DSDV(1,1,i)= (1.-z)*(y-1.);
830 DSDV(1,2,i)= (1.-x)*(z-1.);
831 DSDV(1,3,i)= (1.-x)*(y-1.);
832 DSDV(2,1,i)= (1.-z)*(1.-y);
833 DSDV(2,2,i)= (z-1.)*x;
834 DSDV(2,3,i)= (y-1.)*x;
835 DSDV(3,1,i)= (1.-z)*y;
836 DSDV(3,2,i)= (1.-z)*x;
837 DSDV(3,3,i)=-y*x;
838 DSDV(4,1,i)= y*(z-1.);
839 DSDV(4,2,i)= (1.-z)*(1.-x);
840 DSDV(4,3,i)= y*(x-1.);
841 DSDV(5,1,i)= z*(y-1.);
842 DSDV(5,2,i)= z*(x-1.);
843 DSDV(5,3,i)= (x-1.)*(y-1.);
844 DSDV(6,1,i)= z*(1.-y);
845 DSDV(6,2,i)= -x*z;
846 DSDV(6,3,i)= (1.-y)*x;
847 DSDV(7,1,i)= y*z;
848 DSDV(7,2,i)= x*z;
849 DSDV(7,3,i)= x*y;
850 DSDV(8,1,i)=-y*z;
851 DSDV(8,2,i)= z*(1.-x);
852 DSDV(8,3,i)= y*(1.-x);
853 }
854 #undef NUMSHAPES
855 #undef DIM
856 }
857
858 /************************************************************************************/
859
860 void Finley_Shape_Hex20(int NumV,double* v,double* s,double* dsdv) {
861 #define NUMSHAPES 20
862 #define DIM 3
863 register double x,y,z;
864 int i;
865 #pragma ivdep
866 for (i=0;i<NumV;i++) {
867 x=V(1,i);
868 y=V(2,i);
869 z=V(3,i);
870 S(1,i)=1.+(-3.)*x+(-3.)*y+(-3.)*z+(5.)*x*y+(5.)*x*z+(5.)*y*z+(2.)*x*x+(2.)*y*y+(2.)*z*z+(-2.)*x*x*y+(-2.)*x*x*z+(-2.)*x*y*y+(-2.)*y*y*z+(-2.)*x*z*z+(-2.)*y*z*z+(-7.)*x*y*z+(2.)*x*x*y*z+(2.)*x*y*y*z+(2.)*x*y*z*z;
871 S(2,i)=(-1.)*x+(-1.)*x*y+(-1.)*x*z+(2.)*x*x+(-2.)*x*x*y+(-2.)*x*x*z+(2.)*x*y*y+(2.)*x*z*z+(3.)*x*y*z+(2.)*x*x*y*z+(-2.)*x*y*y*z+(-2.)*x*y*z*z;
872 S(3,i)=(-3.)*x*y+(2.)*x*x*y+(2.)*x*y*y+1.*x*y*z+(-2.)*x*x*y*z+(-2.)*x*y*y*z+(2.)*x*y*z*z;
873 S(4,i)=(-1.)*y+(-1.)*x*y+(-1.)*y*z+(2.)*y*y+(2.)*x*x*y+(-2.)*x*y*y+(-2.)*y*y*z+(2.)*y*z*z+(3.)*x*y*z+(-2.)*x*x*y*z+(2.)*x*y*y*z+(-2.)*x*y*z*z;
874 S(5,i)=(-1.)*z+(-1.)*x*z+(-1.)*y*z+(2.)*z*z+(2.)*x*x*z+(2.)*y*y*z+(-2.)*x*z*z+(-2.)*y*z*z+(3.)*x*y*z+(-2.)*x*x*y*z+(-2.)*x*y*y*z+(2.)*x*y*z*z;
875 S(6,i)=(-3.)*x*z+(2.)*x*x*z+(2.)*x*z*z+1.*x*y*z+(-2.)*x*x*y*z+(2.)*x*y*y*z+(-2.)*x*y*z*z;
876 S(7,i)=(-5.)*x*y*z+(2.)*x*x*y*z+(2.)*x*y*y*z+(2.)*x*y*z*z;
877 S(8,i)=(-3.)*y*z+(2.)*y*y*z+(2.)*y*z*z+1.*x*y*z+(2.)*x*x*y*z+(-2.)*x*y*y*z+(-2.)*x*y*z*z;
878 S(9,i)=(4.)*x+(-4.)*x*y+(-4.)*x*z+(-4.)*x*x+(4.)*x*x*y+(4.)*x*x*z+(4.)*x*y*z+(-4.)*x*x*y*z;
879 S(10,i)=(4.)*x*y+(-4.)*x*y*y+(-4.)*x*y*z+(4.)*x*y*y*z;
880 S(11,i)=(4.)*x*y+(-4.)*x*x*y+(-4.)*x*y*z+(4.)*x*x*y*z;
881 S(12,i)=(4.)*y+(-4.)*x*y+(-4.)*y*z+(-4.)*y*y+(4.)*x*y*y+(4.)*y*y*z+(4.)*x*y*z+(-4.)*x*y*y*z;
882 S(13,i)=(4.)*z+(-4.)*x*z+(-4.)*y*z+(-4.)*z*z+(4.)*x*z*z+(4.)*y*z*z+(4.)*x*y*z+(-4.)*x*y*z*z;
883 S(14,i)=(4.)*x*z+(-4.)*x*z*z+(-4.)*x*y*z+(4.)*x*y*z*z;
884 S(15,i)=(4.)*x*y*z+(-4.)*x*y*z*z;
885 S(16,i)=(4.)*y*z+(-4.)*y*z*z+(-4.)*x*y*z+(4.)*x*y*z*z;
886 S(17,i)=(4.)*x*z+(-4.)*x*x*z+(-4.)*x*y*z+(4.)*x*x*y*z;
887 S(18,i)=(4.)*x*y*z+(-4.)*x*y*y*z;
888 S(19,i)=(4.)*x*y*z+(-4.)*x*x*y*z;
889 S(20,i)=(4.)*y*z+(-4.)*y*y*z+(-4.)*x*y*z+(4.)*x*y*y*z;
890 DSDV(1,1,i)=(-3.)+(5.)*y+(5.)*z+(4.)*x+(-4.)*x*y+(-4.)*x*z+(-2.)*y*y+(-2.)*z*z+(-7.)*y*z+(4.)*x*y*z+(2.)*y*y*z+(2.)*y*z*z;
891 DSDV(2,1,i)=(-1.)+(-1.)*y+(-1.)*z+(4.)*x+(-4.)*x*y+(-4.)*x*z+(2.)*y*y+(2.)*z*z+(3.)*y*z+(4.)*x*y*z+(-2.)*y*y*z+(-2.)*y*z*z;
892 DSDV(3,1,i)=(-3.)*y+(4.)*x*y+(2.)*y*y+1.*y*z+(-4.)*x*y*z+(-2.)*y*y*z+(2.)*y*z*z;
893 DSDV(4,1,i)=(-1.)*y+(4.)*x*y+(-2.)*y*y+(3.)*y*z+(-4.)*x*y*z+(2.)*y*y*z+(-2.)*y*z*z;
894 DSDV(5,1,i)=(-1.)*z+(4.)*x*z+(-2.)*z*z+(3.)*y*z+(-4.)*x*y*z+(-2.)*y*y*z+(2.)*y*z*z;
895 DSDV(6,1,i)=(-3.)*z+(4.)*x*z+(2.)*z*z+1.*y*z+(-4.)*x*y*z+(2.)*y*y*z+(-2.)*y*z*z;
896 DSDV(7,1,i)=(-5.)*y*z+(4.)*x*y*z+(2.)*y*y*z+(2.)*y*z*z;
897 DSDV(8,1,i)=1.*y*z+(4.)*x*y*z+(-2.)*y*y*z+(-2.)*y*z*z;
898 DSDV(9,1,i)=(4.)+(-4.)*y+(-4.)*z+(-8.)*x+(8.)*x*y+(8.)*x*z+(4.)*y*z+(-8.)*x*y*z;
899 DSDV(10,1,i)=(4.)*y+(-4.)*y*y+(-4.)*y*z+(4.)*y*y*z;
900 DSDV(11,1,i)=(4.)*y+(-8.)*x*y+(-4.)*y*z+(8.)*x*y*z;
901 DSDV(12,1,i)=(-4.)*y+(4.)*y*y+(4.)*y*z+(-4.)*y*y*z;
902 DSDV(13,1,i)=(-4.)*z+(4.)*z*z+(4.)*y*z+(-4.)*y*z*z;
903 DSDV(14,1,i)=(4.)*z+(-4.)*z*z+(-4.)*y*z+(4.)*y*z*z;
904 DSDV(15,1,i)=(4.)*y*z+(-4.)*y*z*z;
905 DSDV(16,1,i)=(-4.)*y*z+(4.)*y*z*z;
906 DSDV(17,1,i)=(4.)*z+(-8.)*x*z+(-4.)*y*z+(8.)*x*y*z;
907 DSDV(18,1,i)=(4.)*y*z+(-4.)*y*y*z;
908 DSDV(19,1,i)=(4.)*y*z+(-8.)*x*y*z;
909 DSDV(20,1,i)=(-4.)*y*z+(4.)*y*y*z;
910 DSDV(1,2,i)=(-3.)+(5.)*x+(5.)*z+(4.)*y+(-2.)*x*x+(-4.)*x*y+(-4.)*y*z+(-2.)*z*z+(-7.)*x*z+(2.)*x*x*z+(4.)*x*y*z+(2.)*x*z*z;
911 DSDV(2,2,i)=(-1.)*x+(-2.)*x*x+(4.)*x*y+(3.)*x*z+(2.)*x*x*z+(-4.)*x*y*z+(-2.)*x*z*z;
912 DSDV(3,2,i)=(-3.)*x+(2.)*x*x+(4.)*x*y+1.*x*z+(-2.)*x*x*z+(-4.)*x*y*z+(2.)*x*z*z;
913 DSDV(4,2,i)=(-1.)+(-1.)*x+(-1.)*z+(4.)*y+(2.)*x*x+(-4.)*x*y+(-4.)*y*z+(2.)*z*z+(3.)*x*z+(-2.)*x*x*z+(4.)*x*y*z+(-2.)*x*z*z;
914 DSDV(5,2,i)=(-1.)*z+(4.)*y*z+(-2.)*z*z+(3.)*x*z+(-2.)*x*x*z+(-4.)*x*y*z+(2.)*x*z*z;
915 DSDV(6,2,i)=1.*x*z+(-2.)*x*x*z+(4.)*x*y*z+(-2.)*x*z*z;
916 DSDV(7,2,i)=(-5.)*x*z+(2.)*x*x*z+(4.)*x*y*z+(2.)*x*z*z;
917 DSDV(8,2,i)=(-3.)*z+(4.)*y*z+(2.)*z*z+1.*x*z+(2.)*x*x*z+(-4.)*x*y*z+(-2.)*x*z*z;
918 DSDV(9,2,i)=(-4.)*x+(4.)*x*x+(4.)*x*z+(-4.)*x*x*z;
919 DSDV(10,2,i)=(4.)*x+(-8.)*x*y+(-4.)*x*z+(8.)*x*y*z;
920 DSDV(11,2,i)=(4.)*x+(-4.)*x*x+(-4.)*x*z+(4.)*x*x*z;
921 DSDV(12,2,i)=(4.)+(-4.)*x+(-4.)*z+(-8.)*y+(8.)*x*y+(8.)*y*z+(4.)*x*z+(-8.)*x*y*z;
922 DSDV(13,2,i)=(-4.)*z+(4.)*z*z+(4.)*x*z+(-4.)*x*z*z;
923 DSDV(14,2,i)=(-4.)*x*z+(4.)*x*z*z;
924 DSDV(15,2,i)=(4.)*x*z+(-4.)*x*z*z;
925 DSDV(16,2,i)=(4.)*z+(-4.)*z*z+(-4.)*x*z+(4.)*x*z*z;
926 DSDV(17,2,i)=(-4.)*x*z+(4.)*x*x*z;
927 DSDV(18,2,i)=(4.)*x*z+(-8.)*x*y*z;
928 DSDV(19,2,i)=(4.)*x*z+(-4.)*x*x*z;
929 DSDV(20,2,i)=(4.)*z+(-8.)*y*z+(-4.)*x*z+(8.)*x*y*z;
930 DSDV(1,3,i)=(-3.)+(5.)*x+(5.)*y+(4.)*z+(-2.)*x*x+(-2.)*y*y+(-4.)*x*z+(-4.)*y*z+(-7.)*x*y+(2.)*x*x*y+(2.)*x*y*y+(4.)*x*y*z;
931 DSDV(2,3,i)=(-1.)*x+(-2.)*x*x+(4.)*x*z+(3.)*x*y+(2.)*x*x*y+(-2.)*x*y*y+(-4.)*x*y*z;
932 DSDV(3,3,i)=1.*x*y+(-2.)*x*x*y+(-2.)*x*y*y+(4.)*x*y*z;
933 DSDV(4,3,i)=(-1.)*y+(-2.)*y*y+(4.)*y*z+(3.)*x*y+(-2.)*x*x*y+(2.)*x*y*y+(-4.)*x*y*z;
934 DSDV(5,3,i)=(-1.)+(-1.)*x+(-1.)*y+(4.)*z+(2.)*x*x+(2.)*y*y+(-4.)*x*z+(-4.)*y*z+(3.)*x*y+(-2.)*x*x*y+(-2.)*x*y*y+(4.)*x*y*z;
935 DSDV(6,3,i)=(-3.)*x+(2.)*x*x+(4.)*x*z+1.*x*y+(-2.)*x*x*y+(2.)*x*y*y+(-4.)*x*y*z;
936 DSDV(7,3,i)=(-5.)*x*y+(2.)*x*x*y+(2.)*x*y*y+(4.)*x*y*z;
937 DSDV(8,3,i)=(-3.)*y+(2.)*y*y+(4.)*y*z+1.*x*y+(2.)*x*x*y+(-2.)*x*y*y+(-4.)*x*y*z;
938 DSDV(9,3,i)=(-4.)*x+(4.)*x*x+(4.)*x*y+(-4.)*x*x*y;
939 DSDV(10,3,i)=(-4.)*x*y+(4.)*x*y*y;
940 DSDV(11,3,i)=(-4.)*x*y+(4.)*x*x*y;
941 DSDV(12,3,i)=(-4.)*y+(4.)*y*y+(4.)*x*y+(-4.)*x*y*y;
942 DSDV(13,3,i)=(4.)+(-4.)*x+(-4.)*y+(-8.)*z+(8.)*x*z+(8.)*y*z+(4.)*x*y+(-8.)*x*y*z;
943 DSDV(14,3,i)=(4.)*x+(-8.)*x*z+(-4.)*x*y+(8.)*x*y*z;
944 DSDV(15,3,i)=(4.)*x*y+(-8.)*x*y*z;
945 DSDV(16,3,i)=(4.)*y+(-8.)*y*z+(-4.)*x*y+(8.)*x*y*z;
946 DSDV(17,3,i)=(4.)*x+(-4.)*x*x+(-4.)*x*y+(4.)*x*x*y;
947 DSDV(18,3,i)=(4.)*x*y+(-4.)*x*y*y;
948 DSDV(19,3,i)=(4.)*x*y+(-4.)*x*x*y;
949 DSDV(20,3,i)=(4.)*y+(-4.)*y*y+(-4.)*x*y+(4.)*x*y*y;
950 }
951 #undef NUMSHAPES
952 #undef DIM
953 }
954
955 /************************************************************************************/
956
957 void Finley_Shape_Hex27(int NumV,double* v,double* s,double* dsdv) {
958 #define NUMSHAPES 27
959 #define DIM 3
960 register double x,y,z;
961 int i;
962 #pragma ivdep
963 for (i=0;i<NumV;i++) {
964 x=V(1,i);
965 y=V(2,i);
966 z=V(3,i);
967 S(1,i)= + 1.0 - 3.0*x + 2.0*x*x - 3.0*y + 9.0*x*y - 6.0*x*x*y + 2.0*y*y - 6.0*x*y*y + 4.0*x*x*y*y - 3.0*z + 9.0*x*z - 6.0*x*x*z + 9.0*y*z - 27.0*x*y*z + 18.0*x*x*y*z - 6.0*y*y*z + 18.0*x*y*y*z - 12.0*x*x*y*y*z + 2.0*z*z - 6.0*x*z*z + 4.0*x*x*z*z - 6.0*y*z*z + 18.0*x*y*z*z - 12.0*x*x*y*z*z + 4.0*y*y*z*z - 12.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
968 S(2,i)= - 1.0*x + 2.0*x*x + 3.0*x*y - 6.0*x*x*y - 2.0*x*y*y + 4.0*x*x*y*y + 3.0*x*z - 6.0*x*x*z - 9.0*x*y*z + 18.0*x*x*y*z + 6.0*x*y*y*z - 12.0*x*x*y*y*z - 2.0*x*z*z + 4.0*x*x*z*z + 6.0*x*y*z*z - 12.0*x*x*y*z*z - 4.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
969 S(3,i)= + 1.0*x*y - 2.0*x*x*y - 2.0*x*y*y + 4.0*x*x*y*y - 3.0*x*y*z + 6.0*x*x*y*z + 6.0*x*y*y*z - 12.0*x*x*y*y*z + 2.0*x*y*z*z - 4.0*x*x*y*z*z - 4.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
970 S(4,i)= - 1.0*y + 3.0*x*y - 2.0*x*x*y + 2.0*y*y - 6.0*x*y*y + 4.0*x*x*y*y + 3.0*y*z - 9.0*x*y*z + 6.0*x*x*y*z - 6.0*y*y*z + 18.0*x*y*y*z - 12.0*x*x*y*y*z - 2.0*y*z*z + 6.0*x*y*z*z - 4.0*x*x*y*z*z + 4.0*y*y*z*z - 12.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
971 S(5,i)= - 1.0*z + 3.0*x*z - 2.0*x*x*z + 3.0*y*z - 9.0*x*y*z + 6.0*x*x*y*z - 2.0*y*y*z + 6.0*x*y*y*z - 4.0*x*x*y*y*z + 2.0*z*z - 6.0*x*z*z + 4.0*x*x*z*z - 6.0*y*z*z + 18.0*x*y*z*z - 12.0*x*x*y*z*z + 4.0*y*y*z*z - 12.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
972 S(6,i)= + 1.0*x*z - 2.0*x*x*z - 3.0*x*y*z + 6.0*x*x*y*z + 2.0*x*y*y*z - 4.0*x*x*y*y*z - 2.0*x*z*z + 4.0*x*x*z*z + 6.0*x*y*z*z - 12.0*x*x*y*z*z - 4.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
973 S(7,i)= - 1.0*x*y*z + 2.0*x*x*y*z + 2.0*x*y*y*z - 4.0*x*x*y*y*z + 2.0*x*y*z*z - 4.0*x*x*y*z*z - 4.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
974 S(8,i)= + 1.0*y*z - 3.0*x*y*z + 2.0*x*x*y*z - 2.0*y*y*z + 6.0*x*y*y*z - 4.0*x*x*y*y*z - 2.0*y*z*z + 6.0*x*y*z*z - 4.0*x*x*y*z*z + 4.0*y*y*z*z - 12.0*x*y*y*z*z + 8.0*x*x*y*y*z*z;
975 S(9,i)= + 4.0*x - 4.0*x*x - 12.0*x*y + 12.0*x*x*y + 8.0*x*y*y - 8.0*x*x*y*y - 12.0*x*z + 12.0*x*x*z + 36.0*x*y*z - 36.0*x*x*y*z - 24.0*x*y*y*z + 24.0*x*x*y*y*z + 8.0*x*z*z - 8.0*x*x*z*z - 24.0*x*y*z*z + 24.0*x*x*y*z*z + 16.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
976 S(10,i)= - 4.0*x*y + 8.0*x*x*y + 4.0*x*y*y - 8.0*x*x*y*y + 12.0*x*y*z - 24.0*x*x*y*z - 12.0*x*y*y*z + 24.0*x*x*y*y*z - 8.0*x*y*z*z + 16.0*x*x*y*z*z + 8.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
977 S(11,i)= - 4.0*x*y + 4.0*x*x*y + 8.0*x*y*y - 8.0*x*x*y*y + 12.0*x*y*z - 12.0*x*x*y*z - 24.0*x*y*y*z + 24.0*x*x*y*y*z - 8.0*x*y*z*z + 8.0*x*x*y*z*z + 16.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
978 S(12,i)= + 4.0*y - 12.0*x*y + 8.0*x*x*y - 4.0*y*y + 12.0*x*y*y - 8.0*x*x*y*y - 12.0*y*z + 36.0*x*y*z - 24.0*x*x*y*z + 12.0*y*y*z - 36.0*x*y*y*z + 24.0*x*x*y*y*z + 8.0*y*z*z - 24.0*x*y*z*z + 16.0*x*x*y*z*z - 8.0*y*y*z*z + 24.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
979 S(13,i)= + 4.0*z - 12.0*x*z + 8.0*x*x*z - 12.0*y*z + 36.0*x*y*z - 24.0*x*x*y*z + 8.0*y*y*z - 24.0*x*y*y*z + 16.0*x*x*y*y*z - 4.0*z*z + 12.0*x*z*z - 8.0*x*x*z*z + 12.0*y*z*z - 36.0*x*y*z*z + 24.0*x*x*y*z*z - 8.0*y*y*z*z + 24.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
980 S(14,i)= - 4.0*x*z + 8.0*x*x*z + 12.0*x*y*z - 24.0*x*x*y*z - 8.0*x*y*y*z + 16.0*x*x*y*y*z + 4.0*x*z*z - 8.0*x*x*z*z - 12.0*x*y*z*z + 24.0*x*x*y*z*z + 8.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
981 S(15,i)= + 4.0*x*y*z - 8.0*x*x*y*z - 8.0*x*y*y*z + 16.0*x*x*y*y*z - 4.0*x*y*z*z + 8.0*x*x*y*z*z + 8.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
982 S(16,i)= - 4.0*y*z + 12.0*x*y*z - 8.0*x*x*y*z + 8.0*y*y*z - 24.0*x*y*y*z + 16.0*x*x*y*y*z + 4.0*y*z*z - 12.0*x*y*z*z + 8.0*x*x*y*z*z - 8.0*y*y*z*z + 24.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
983 S(17,i)= - 4.0*x*z + 4.0*x*x*z + 12.0*x*y*z - 12.0*x*x*y*z - 8.0*x*y*y*z + 8.0*x*x*y*y*z + 8.0*x*z*z - 8.0*x*x*z*z - 24.0*x*y*z*z + 24.0*x*x*y*z*z + 16.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
984 S(18,i)= + 4.0*x*y*z - 8.0*x*x*y*z - 4.0*x*y*y*z + 8.0*x*x*y*y*z - 8.0*x*y*z*z + 16.0*x*x*y*z*z + 8.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
985 S(19,i)= + 4.0*x*y*z - 4.0*x*x*y*z - 8.0*x*y*y*z + 8.0*x*x*y*y*z - 8.0*x*y*z*z + 8.0*x*x*y*z*z + 16.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
986 S(20,i)= - 4.0*y*z + 12.0*x*y*z - 8.0*x*x*y*z + 4.0*y*y*z - 12.0*x*y*y*z + 8.0*x*x*y*y*z + 8.0*y*z*z - 24.0*x*y*z*z + 16.0*x*x*y*z*z - 8.0*y*y*z*z + 24.0*x*y*y*z*z - 16.0*x*x*y*y*z*z;
987 S(21,i)= + 16.0*x*y - 16.0*x*x*y - 16.0*x*y*y + 16.0*x*x*y*y - 48.0*x*y*z + 48.0*x*x*y*z + 48.0*x*y*y*z - 48.0*x*x*y*y*z + 32.0*x*y*z*z - 32.0*x*x*y*z*z - 32.0*x*y*y*z*z + 32.0*x*x*y*y*z*z;
988 S(22,i)= + 16.0*x*z - 16.0*x*x*z - 48.0*x*y*z + 48.0*x*x*y*z + 32.0*x*y*y*z - 32.0*x*x*y*y*z - 16.0*x*z*z + 16.0*x*x*z*z + 48.0*x*y*z*z - 48.0*x*x*y*z*z - 32.0*x*y*y*z*z + 32.0*x*x*y*y*z*z;
989 S(23,i)= - 16.0*x*y*z + 32.0*x*x*y*z + 16.0*x*y*y*z - 32.0*x*x*y*y*z + 16.0*x*y*z*z - 32.0*x*x*y*z*z - 16.0*x*y*y*z*z + 32.0*x*x*y*y*z*z;
990 S(24,i)= - 16.0*x*y*z + 16.0*x*x*y*z + 32.0*x*y*y*z - 32.0*x*x*y*y*z + 16.0*x*y*z*z - 16.0*x*x*y*z*z - 32.0*x*y*y*z*z + 32.0*x*x*y*y*z*z;
991 S(25,i)= + 16.0*y*z - 48.0*x*y*z + 32.0*x*x*y*z - 16.0*y*y*z + 48.0*x*y*y*z - 32.0*x*x*y*y*z - 16.0*y*z*z + 48.0*x*y*z*z - 32.0*x*x*y*z*z + 16.0*y*y*z*z - 48.0*x*y*y*z*z + 32.0*x*x*y*y*z*z;
992 S(26,i)= - 16.0*x*y*z + 16.0*x*x*y*z + 16.0*x*y*y*z - 16.0*x*x*y*y*z + 32.0*x*y*z*z - 32.0*x*x*y*z*z - 32.0*x*y*y*z*z + 32.0*x*x*y*y*z*z;
993 S(27,i)= + 64.0*x*y*z - 64.0*x*x*y*z - 64.0*x*y*y*z + 64.0*x*x*y*y*z - 64.0*x*y*z*z + 64.0*x*x*y*z*z + 64.0*x*y*y*z*z - 64.0*x*x*y*y*z*z;
994 DSDV(1,1,i)= - 3.0 + 4.0*x + 9.0*y - 12.0*x*y - 6.0*y*y + 8.0*x*y*y + 9.0*z - 12.0*x*z - 27.0*y*z + 36.0*x*y*z + 18.0*y*y*z - 24.0*x*y*y*z - 6.0*z*z + 8.0*x*z*z + 18.0*y*z*z - 24.0*x*y*z*z - 12.0*y*y*z*z + 16.0*x*y*y*z*z;
995 DSDV(1,2,i)= - 3.0 + 9.0*x - 6.0*x*x + 4.0*y - 12.0*x*y + 8.0*x*x*y + 9.0*z - 27.0*x*z + 18.0*x*x*z - 12.0*y*z + 36.0*x*y*z - 24.0*x*x*y*z - 6.0*z*z + 18.0*x*z*z - 12.0*x*x*z*z + 8.0*y*z*z - 24.0*x*y*z*z + 16.0*x*x*y*z*z;
996 DSDV(1,3,i)= - 3.0 + 9.0*x - 6.0*x*x + 9.0*y - 27.0*x*y + 18.0*x*x*y - 6.0*y*y + 18.0*x*y*y - 12.0*x*x*y*y + 4.0*z - 12.0*x*z + 8.0*x*x*z - 12.0*y*z + 36.0*x*y*z - 24.0*x*x*y*z + 8.0*y*y*z - 24.0*x