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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3171 - (show annotations)
Fri Sep 10 00:31:11 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 20970 byte(s)
Remove some contact jacobean functions.
Assemble_jacobeans_2D_M1D_E2D finished first pass

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 #include "Assemble.h"
16 #include "Util.h"
17 #ifdef _OPENMP
18 #include <omp.h>
19 #endif
20
21
22 /* input:
23
24 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 /**************************************************************/
46 /* */
47 /* Jacobean 2D with area element */
48 /* */
49 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 #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 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 {
100 #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 }
116 volume[INDEX2(q,e,3)]=ABS(D)*quadweight;
117 }
118 }
119 }
120 }
121 }
122 #undef DIM
123 #undef LOCDIM
124 #undef DTDXSET
125 #undef COMPDXDV0
126 #undef COMPDXDV1
127 }
128 /**************************************************************/
129 /* */
130 /* Jacobean 1D manifold in 2D and 1D elements */
131 /* */
132 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,/*double* bogus_QuadWeights,
133 dim_t bogus_numShape,*/ dim_t numElements, dim_t numNodes, index_t* nodes,
134 /*double* bogus_DSDv, dim_t bogus_numTest, double* DTDv,*/
135 double* dTdX, double* volume, index_t* element_id) {
136 #define DIM 2
137 #define LOCDIM 1
138 register int e;
139 char error_msg[LenErrorMsg_MAX];
140 const dim_t numTest=2;
141 double quadweight=(numQuad==1)?1.0:0.5;
142 // numQuad is 1 or 2
143 #pragma omp parallel
144 {
145 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
146 #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
147 for(e=0;e<numElements;e++){
148 dXdv00=0;
149 dXdv10=0;
150 dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
151 dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
152 D=dXdv00*dXdv00+dXdv10*dXdv10;
153 if (D==0.) {
154 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
155 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
156 } else {
157 invD=1./D;
158 dvdX00=dXdv00*invD;
159 dvdX01=dXdv10*invD;
160 // The number of quad points is 1 or 2
161
162
163 dTdX[INDEX4(0,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
164 dTdX[INDEX4(0,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
165 dTdX[INDEX4(1,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
166 dTdX[INDEX4(1,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
167 volume[INDEX2(0,e,numQuad)]=sqrt(D)*quadweight;
168
169 if (numQuad==2)
170 {
171
172 dTdX[INDEX4(0,0,1,e,numTest,DIM,numQuad)]=dvdX00;
173 dTdX[INDEX4(0,1,1,e,numTest,DIM,numQuad)]=dvdX01;
174 dTdX[INDEX4(1,0,1,e,numTest,DIM,numQuad)]=dvdX00;
175 dTdX[INDEX4(1,1,1,e,numTest,DIM,numQuad)]=dvdX01;
176 volume[INDEX2(1,e,numQuad)]=sqrt(D)*quadweight;
177
178 }
179 }
180 }
181
182 }
183 #undef DIM
184 #undef LOCDIM
185 }
186
187 /**************************************************************/
188 /* */
189 /* Jacobean 1D manifold in 2D and 2D elements */
190 /* */
191 void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
192 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
193 double* DSDv, dim_t numTest,double* DTDv,
194 double* dTdX, double* volume, index_t* element_id) {
195 #define DIM 2
196 #define LOCDIM 2
197 register int e,q,s;
198 char error_msg[LenErrorMsg_MAX];
199 #pragma omp parallel
200 {
201 register double dXdv00,dXdv10,dXdv01,dXdv11,
202 dvdX00,dvdX10,dvdX01,dvdX11, D,invD,
203 X0_loc, X1_loc;
204 #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
205 for(e=0;e<numElements;e++){
206 for (q=0;q<numQuad;q++) {
207 dXdv00=0;
208 dXdv10=0;
209 dXdv01=0;
210 dXdv11=0;
211 for (s=0;s<numShape; s++) {
212 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
213 X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
214 dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
215 dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
216 dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
217 dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
218 }
219 D = dXdv00*dXdv11 - dXdv01*dXdv10;
220 if (D==0.) {
221 sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
222 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
223 } else {
224 invD=1./D;
225 dvdX00= dXdv11*invD;
226 dvdX10=-dXdv10*invD;
227 dvdX01=-dXdv01*invD;
228 dvdX11= dXdv00*invD;
229
230 for (s=0;s<numTest; s++) {
231 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;
232 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;
233 }
234 }
235 volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
236 }
237 }
238
239 }
240 #undef DIM
241 #undef LOCDIM
242 }
243
244 /**************************************************************/
245 /* */
246 /* Jacobean 3D */
247 /* */
248 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
249 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
250 double* DSDv, dim_t numTest, double* DTDv,
251 double* dTdX, double* volume, index_t* element_id) {
252 #define DIM 3
253 #define LOCDIM 3
254 int e,q,s;
255 char error_msg[LenErrorMsg_MAX];
256 #pragma omp parallel
257 {
258 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
259 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
260 X0_loc,X1_loc,X2_loc;
261
262 #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)
263 for(e=0;e<numElements;e++){
264 for (q=0;q<numQuad;q++) {
265 dXdv00=0;
266 dXdv10=0;
267 dXdv20=0;
268 dXdv01=0;
269 dXdv11=0;
270 dXdv21=0;
271 dXdv02=0;
272 dXdv12=0;
273 dXdv22=0;
274 for (s=0;s<numShape; s++) {
275 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
276 X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
277 X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
278 dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
279 dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
280 dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
281 dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
282 dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
283 dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
284 dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
285 dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
286 dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
287 }
288 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
289 if (D==0.) {
290 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
291 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
292 } else {
293 invD=1./D;
294 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
295 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
296 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
297 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
298 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
299 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
300 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
301 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
302 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
303
304 for (s=0;s<numTest; s++) {
305 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+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
306 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+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
307 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+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
308 }
309 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
310 }
311 }
312 }
313
314 }
315 #undef DIM
316 #undef LOCDIM
317 }
318 /**************************************************************/
319 /* */
320 /* Jacobean 2D manifold in 3D with 3D elements */
321 /* */
322 void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
323 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
324 double* DSDv, dim_t numTest,double* DTDv,
325 double* dTdX, double* volume, index_t* element_id) {
326 #define DIM 3
327 #define LOCDIM 3
328 register int e,q,s;
329 char error_msg[LenErrorMsg_MAX];
330 #pragma omp parallel
331 {
332 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,
333 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
334 X0_loc, X1_loc, X2_loc;
335 #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)
336 for(e=0;e<numElements;e++){
337 for (q=0;q<numQuad;q++) {
338 dXdv00=0;
339 dXdv10=0;
340 dXdv20=0;
341 dXdv01=0;
342 dXdv11=0;
343 dXdv21=0;
344 dXdv02=0;
345 dXdv12=0;
346 dXdv22=0;
347 for (s=0;s<numShape; s++) {
348 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
349 X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
350 X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
351 dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
352 dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
353 dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
354 dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
355 dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
356 dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
357 dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
358 dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
359 dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
360 }
361 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
362 if (D==0.) {
363 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
364 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
365 } else {
366 invD=1./D;
367 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
368 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
369 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
370 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
371 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
372 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
373 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
374 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
375 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
376
377 for (s=0;s<numTest; s++) {
378 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=
379 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
380 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=
381 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
382 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=
383 DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
384 }
385 }
386 m0=dXdv10*dXdv21-dXdv20*dXdv11;
387 m1=dXdv20*dXdv01-dXdv00*dXdv21;
388 m2=dXdv00*dXdv11-dXdv10*dXdv01;
389 volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
390 }
391 }
392
393 }
394 #undef DIM
395 #undef LOCDIM
396 }
397 /**************************************************************/
398 /* */
399 /* Jacobean 2D manifold in 3D with 2D elements */
400 /* */
401 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
402 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
403 double* DSDv, dim_t numTest,double* DTDv,
404 double* dTdX, double* volume, index_t* element_id) {
405 #define DIM 3
406 #define LOCDIM 2
407 register int e,q,s;
408 char error_msg[LenErrorMsg_MAX];
409 #pragma omp parallel
410 {
411 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
412 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
413 X0_loc, X1_loc, X2_loc;
414 #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)
415 for(e=0;e<numElements;e++){
416 for (q=0;q<numQuad;q++) {
417 dXdv00=0;
418 dXdv10=0;
419 dXdv20=0;
420 dXdv01=0;
421 dXdv11=0;
422 dXdv21=0;
423 for (s=0;s<numShape; s++) {
424 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
425 X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
426 X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
427 dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
428 dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
429 dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
430 dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
431 dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
432 dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
433 }
434 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
435 m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
436 m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
437 D=m00*m11-m01*m01;
438 if (D==0.) {
439 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
440 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
441 } else {
442 invD=1./D;
443 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
444 dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
445 dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
446 dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
447 dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
448 dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
449 for (s=0;s<numTest; s++) {
450 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;
451 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;
452 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;
453 }
454 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
455 }
456 }
457 }
458
459 }
460 #undef DIM
461 #undef LOCDIM
462 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26