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

Annotation of /branches/domexper/dudley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3172 - (hide annotations)
Fri Sep 10 01:38:04 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 17749 byte(s)
Jacobeans_3D done first pass

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 jgs 82
14 ksteube 1811
15 jgs 82 #include "Assemble.h"
16     #include "Util.h"
17     #ifdef _OPENMP
18     #include <omp.h>
19     #endif
20    
21 gross 777
22 gross 2748 /* input:
23 gross 777
24 gross 2748 double* coordinates[DIM*(*)]
25     dim_t numQuad
26     double* QuadWeights[numQuad]
27     dim_t numShape
28     dim_t numElements
29     dim_t numNodes
30     index_t* nodes[numNodes*numElements] where NUMSIDES*numShape<=numNodes
31     double* DSDv[numShape*DIM*numQuad]
32     dim_t numTest
33     double* DTDv[LOCDIM*numTest*numQuad]
34     index_t* element_id[numElements]
35    
36     output:
37    
38     double* dTdX[DIM*numTest*NUMSIDES*numQuad*numElements]
39     double* volume[numQuad*numElements]
40    
41     */
42    
43     #define SCALING(_nsub_,_dim_) pow(1./(double)(_nsub_),1./(double)(_dim_))
44    
45 jgs 82 /**************************************************************/
46 gross 776 /* */
47 gross 777 /* Jacobean 2D with area element */
48 gross 776 /* */
49 jfenwick 3170 void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
50     double* dTdX, double* volume, index_t* element_id)
51     {
52     #define DIM 2
53     #define LOCDIM 2
54     register int e,q,s;
55     char error_msg[LenErrorMsg_MAX];
56     const double quadweight=(numQuad==1)?1./2:1./6; /* numQuad is 1 or 3 */
57     const dim_t numTest=3; // hoping this is used in constant folding
58     static const int DTDV[3][2]={{-1,-1}, {1,0}, {0,1}}; // if constant folding is not applied here will need to write
59     // out in full
60     #pragma omp parallel
61     {
62     register double dXdv00,dXdv10,dXdv01,dXdv11,
63     dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
64     #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
65     for(e=0;e<numElements;e++)
66     {
67 jfenwick 3169 #define COMPDXDV0(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
68     coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
69     coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
70    
71     #define COMPDXDV1(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
72     coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
73     coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
74    
75 jfenwick 3170 dXdv00=0;
76     dXdv10=0;
77     dXdv01=0;
78     dXdv11=0;
79     dXdv00=COMPDXDV0(0);
80     dXdv10=COMPDXDV0(1);
81     dXdv01=COMPDXDV1(0);
82     dXdv11=COMPDXDV1(1);
83     D = dXdv00*dXdv11 - dXdv01*dXdv10;
84     if (D==0.)
85     {
86     sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
87     Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
88     }
89     else
90     {
91     invD=1./D;
92     dvdX00= dXdv11*invD;
93     dvdX10=-dXdv10*invD;
94     dvdX01=-dXdv01*invD;
95     dvdX11= dXdv00*invD;
96     if (numQuad==1)
97     {
98     for (s=0;s<3; s++)
99 jfenwick 3169 {
100 jfenwick 3170 #define DTDXSET(P,Q) dTdX[INDEX4(s,P,Q,e,numTest,DIM,numQuad)] = DTDV[s][0]*dvdX0##P+DTDV[s][1]*dvdX1##P
101    
102     DTDXSET(0,0);
103     DTDXSET(1,0);
104     }
105     volume[INDEX2(0,e,1)]=ABS(D)*quadweight;
106     }
107     else // numQuad==3
108     {
109     for (q=0;q<numTest;++q) // relying on unroll loops to optimise this
110     {
111     for (s=0;s<3; s++)
112     {
113     DTDXSET(0,q);
114     DTDXSET(1,q);
115 jfenwick 3169 }
116 jfenwick 3170 volume[INDEX2(q,e,3)]=ABS(D)*quadweight;
117 jfenwick 3169 }
118     }
119 gross 776 }
120 jfenwick 3169 }
121 gross 776 }
122     #undef DIM
123 gross 781 #undef LOCDIM
124 jfenwick 3170 #undef DTDXSET
125 jfenwick 3169 #undef COMPDXDV0
126     #undef COMPDXDV1
127 gross 776 }
128     /**************************************************************/
129     /* */
130 gross 777 /* Jacobean 1D manifold in 2D and 1D elements */
131 gross 776 /* */
132 jfenwick 3172 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,
133     dim_t numElements, dim_t numNodes, index_t* nodes,
134 gross 777 double* dTdX, double* volume, index_t* element_id) {
135     #define DIM 2
136 gross 781 #define LOCDIM 1
137 jfenwick 3171 register int e;
138 gross 776 char error_msg[LenErrorMsg_MAX];
139 jfenwick 3171 const dim_t numTest=2;
140     double quadweight=(numQuad==1)?1.0:0.5;
141     // numQuad is 1 or 2
142 gross 776 #pragma omp parallel
143     {
144 jfenwick 3171 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
145 caltinay 2750 #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
146 gross 776 for(e=0;e<numElements;e++){
147     dXdv00=0;
148     dXdv10=0;
149 jfenwick 3171 dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
150     dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
151 gross 777 D=dXdv00*dXdv00+dXdv10*dXdv10;
152 jfenwick 3171 if (D==0.) {
153 gross 777 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
154 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
155 jfenwick 3171 } else {
156 gross 777 invD=1./D;
157     dvdX00=dXdv00*invD;
158     dvdX01=dXdv10*invD;
159 jfenwick 3171 // The number of quad points is 1 or 2
160 gross 777
161 jfenwick 3170
162 jfenwick 3171 dTdX[INDEX4(0,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
163     dTdX[INDEX4(0,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
164     dTdX[INDEX4(1,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
165     dTdX[INDEX4(1,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
166     volume[INDEX2(0,e,numQuad)]=sqrt(D)*quadweight;
167    
168     if (numQuad==2)
169     {
170    
171     dTdX[INDEX4(0,0,1,e,numTest,DIM,numQuad)]=dvdX00;
172     dTdX[INDEX4(0,1,1,e,numTest,DIM,numQuad)]=dvdX01;
173     dTdX[INDEX4(1,0,1,e,numTest,DIM,numQuad)]=dvdX00;
174     dTdX[INDEX4(1,1,1,e,numTest,DIM,numQuad)]=dvdX01;
175     volume[INDEX2(1,e,numQuad)]=sqrt(D)*quadweight;
176    
177     }
178 gross 781 }
179     }
180    
181     }
182     #undef DIM
183     #undef LOCDIM
184     }
185 jfenwick 3171
186 gross 781 /**************************************************************/
187     /* */
188 gross 777 /* Jacobean 3D */
189 gross 776 /* */
190 jfenwick 3172 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
191 gross 777 double* dTdX, double* volume, index_t* element_id) {
192 gross 776 #define DIM 3
193 gross 781 #define LOCDIM 3
194 gross 853 int e,q,s;
195 gross 776 char error_msg[LenErrorMsg_MAX];
196 jfenwick 3172 // numQuad is 1 or 4
197     const dim_t numShape=4, numTest=4;
198     const double quadweight=(numQuad==1)?1./6:1./24;
199     const double DTDV[4][3]={{-1, -1, -1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
200 gross 776 #pragma omp parallel
201     {
202 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
203 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
204     X0_loc,X1_loc,X2_loc;
205    
206 caltinay 2750 #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22,D,invD,X0_loc,X1_loc,X2_loc) schedule(static)
207 gross 776 for(e=0;e<numElements;e++){
208     dXdv00=0;
209     dXdv10=0;
210     dXdv20=0;
211     dXdv01=0;
212     dXdv11=0;
213     dXdv21=0;
214 gross 777 dXdv02=0;
215     dXdv12=0;
216     dXdv22=0;
217 gross 776 for (s=0;s<numShape; s++) {
218 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
219     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
220     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
221 jfenwick 3172 dXdv00+=X0_loc*DTDV[s][0];
222     dXdv10+=X1_loc*DTDV[s][0];
223     dXdv20+=X2_loc*DTDV[s][0];
224     dXdv01+=X0_loc*DTDV[s][1];
225     dXdv11+=X1_loc*DTDV[s][1];
226     dXdv21+=X2_loc*DTDV[s][1];
227     dXdv02+=X0_loc*DTDV[s][2];
228     dXdv12+=X1_loc*DTDV[s][2];
229     dXdv22+=X2_loc*DTDV[s][2];
230 gross 776 }
231 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
232 gross 776 if (D==0.) {
233 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
234 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
235 gross 776 } else {
236     invD=1./D;
237 gross 777 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
238     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
239     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
240     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
241     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
242     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
243     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
244     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
245     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
246    
247 jfenwick 3172 for (q=0;q<numQuad;q++) {
248    
249    
250 gross 776 for (s=0;s<numTest; s++) {
251 jfenwick 3172 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX00+DTDV[s][1]*dvdX10+DTDV[s][2]*dvdX20;
252     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX01+DTDV[s][1]*dvdX11+DTDV[s][2]*dvdX21;
253     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX02+DTDV[s][1]*dvdX12+DTDV[s][2]*dvdX22;
254 gross 776 }
255 jfenwick 3172 volume[INDEX2(q,e,numQuad)]=ABS(D)*quadweight;
256 gross 776 }
257     }
258     }
259 jgs 82
260 gross 776 }
261     #undef DIM
262 gross 781 #undef LOCDIM
263 gross 776 }
264     /**************************************************************/
265     /* */
266 gross 777 /* Jacobean 2D manifold in 3D with 3D elements */
267 gross 776 /* */
268 gross 777 void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
269     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
270     double* DSDv, dim_t numTest,double* DTDv,
271     double* dTdX, double* volume, index_t* element_id) {
272 gross 776 #define DIM 3
273 gross 781 #define LOCDIM 3
274 gross 776 register int e,q,s;
275     char error_msg[LenErrorMsg_MAX];
276     #pragma omp parallel
277     {
278 gross 781 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,
279 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
280     X0_loc, X1_loc, X2_loc;
281 caltinay 2750 #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,X0_loc, X1_loc, X2_loc) schedule(static)
282 gross 776 for(e=0;e<numElements;e++){
283     for (q=0;q<numQuad;q++) {
284     dXdv00=0;
285     dXdv10=0;
286     dXdv20=0;
287 gross 777 dXdv01=0;
288     dXdv11=0;
289     dXdv21=0;
290     dXdv02=0;
291     dXdv12=0;
292     dXdv22=0;
293 gross 776 for (s=0;s<numShape; s++) {
294 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
295     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
296     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
297     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
298     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
299     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
300     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
301     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
302     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
303     dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
304     dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
305     dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
306 gross 776 }
307 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
308 gross 776 if (D==0.) {
309 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
310 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
311 gross 776 } else {
312     invD=1./D;
313 gross 777 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
314     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
315     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
316     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
317     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
318     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
319     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
320     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
321     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
322    
323 gross 776 for (s=0;s<numTest; s++) {
324 gross 781 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=
325     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
326     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=
327     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
328     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=
329     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
330 gross 776 }
331     }
332 gross 777 m0=dXdv10*dXdv21-dXdv20*dXdv11;
333     m1=dXdv20*dXdv01-dXdv00*dXdv21;
334     m2=dXdv00*dXdv11-dXdv10*dXdv01;
335     volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
336 gross 776 }
337 jgs 82 }
338    
339 gross 776 }
340     #undef DIM
341 gross 781 #undef LOCDIM
342 gross 776 }
343     /**************************************************************/
344     /* */
345 gross 777 /* Jacobean 2D manifold in 3D with 2D elements */
346 gross 776 /* */
347 gross 777 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
348     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
349     double* DSDv, dim_t numTest,double* DTDv,
350     double* dTdX, double* volume, index_t* element_id) {
351     #define DIM 3
352 gross 781 #define LOCDIM 2
353 gross 776 register int e,q,s;
354     char error_msg[LenErrorMsg_MAX];
355     #pragma omp parallel
356 jgs 82 {
357 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
358 gross 853 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
359     X0_loc, X1_loc, X2_loc;
360 caltinay 2750 #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD, X0_loc, X1_loc, X2_loc) schedule(static)
361 gross 776 for(e=0;e<numElements;e++){
362     for (q=0;q<numQuad;q++) {
363     dXdv00=0;
364     dXdv10=0;
365 gross 777 dXdv20=0;
366     dXdv01=0;
367     dXdv11=0;
368     dXdv21=0;
369 gross 776 for (s=0;s<numShape; s++) {
370 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
371     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
372     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
373     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
374     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
375     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
376     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
377     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
378     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
379 gross 776 }
380 gross 777 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
381     m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
382     m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
383     D=m00*m11-m01*m01;
384 gross 776 if (D==0.) {
385 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
386 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
387 gross 776 } else {
388     invD=1./D;
389 gross 777 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
390     dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
391     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
392     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
393     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
394     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
395 gross 776 for (s=0;s<numTest; s++) {
396 gross 781 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10;
397     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11;
398     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12;
399 gross 776 }
400     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
401     }
402     }
403     }
404 jgs 82
405     }
406 gross 776 #undef DIM
407 gross 781 #undef LOCDIM
408 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26