/[escript]/temp/finley/src/Assemble_jacobeans.c
ViewVC logotype

Contents of /temp/finley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 777 - (show annotations)
Wed Jul 12 08:54:45 2006 UTC (13 years, 3 months ago) by gross
Original Path: trunk/finley/src/Assemble_jacobeans.c
File MIME type: text/plain
File size: 21708 byte(s)
integration with persistent jacobeans is running now
1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13
14 /**************************************************************/
15
16 /* Author: gross@access.edu.au */
17 /* Version: $Id$ */
18
19 /**************************************************************/
20
21 #include "Assemble.h"
22 #include "Util.h"
23 #ifdef _OPENMP
24 #include <omp.h>
25 #endif
26
27
28
29 /**************************************************************/
30 /* */
31 /* Jacobean 1D */
32 /* */
33 void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
34 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
35 double* DSDv, dim_t numTest,double* DTDv,
36 double* dTdX, double* volume, index_t* element_id) {
37 #define DIM 1
38 #define LOCALDIM 1
39 register int e,q,s;
40 char error_msg[LenErrorMsg_MAX];
41 #pragma omp parallel
42 {
43 register double D,invD;
44 double X0[numShape];
45 #pragma omp for private(e,q,s) schedule(static)
46 for(e=0;e<numElements;e++){
47 for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
48 for (q=0;q<numQuad;q++) {
49 D=0;
50 for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
51 if (D==0.) {
52 sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);
53 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
54 } else {
55 invD=1./D;
56 for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;
57 }
58 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
59 }
60 }
61
62 }
63 #undef DIM
64 #undef LOCALDIM
65
66 }
67 /**************************************************************/
68 /* */
69 /* Jacobean 2D with area element */
70 /* */
71 void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
72 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
73 double* DSDv, dim_t numTest,double* DTDv,
74 double* dTdX, double* volume, index_t* element_id) {
75 #define DIM 2
76 #define LOCALDIM 2
77 register int e,q,s;
78 char error_msg[LenErrorMsg_MAX];
79 #pragma omp parallel
80 {
81 register double dXdv00,dXdv10,dXdv01,dXdv11,
82 dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
83 double X0[numShape], X1[numShape];
84 #pragma omp for private(e,q,s) schedule(static)
85 for(e=0;e<numElements;e++){
86 for (s=0;s<numShape; s++) {
87 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
88 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
89 }
90 for (q=0;q<numQuad;q++) {
91 dXdv00=0;
92 dXdv10=0;
93 dXdv01=0;
94 dXdv11=0;
95 for (s=0;s<numShape; s++) {
96 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
97 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
98 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
99 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
100 }
101 D = dXdv00*dXdv11 - dXdv01*dXdv10;
102 if (D==0.) {
103 sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
104 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
105 } else {
106 invD=1./D;
107 dvdX00= dXdv11*invD;
108 dvdX10=-dXdv10*invD;
109 dvdX01=-dXdv01*invD;
110 dvdX11= dXdv00*invD;
111
112 for (s=0;s<numTest; s++) {
113 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
114 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
115 }
116 }
117 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
118 }
119 }
120
121 }
122 #undef DIM
123 #undef LOCALDIM
124 }
125 /**************************************************************/
126 /* */
127 /* Jacobean 1D manifold in 2D and 1D elements */
128 /* */
129 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,double* QuadWeights,
130 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
131 double* DSDv, dim_t numTest,double* DTDv,
132 double* dTdX, double* volume, index_t* element_id) {
133 #define DIM 2
134 #define LOCALDIM 1
135 register int e,q,s;
136 char error_msg[LenErrorMsg_MAX];
137 #pragma omp parallel
138 {
139 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
140 double X0[numShape], X1[numShape];
141 #pragma omp for private(e,q,s) schedule(static)
142 for(e=0;e<numElements;e++){
143 for (s=0;s<numShape; s++) {
144 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
145 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
146 }
147 for (q=0;q<numQuad;q++) {
148 dXdv00=0;
149 dXdv10=0;
150 for (s=0;s<numShape; s++) {
151 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
152 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
153 }
154 D=dXdv00*dXdv00+dXdv10*dXdv10;
155 if (D==0.) {
156 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
157 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
158 } else {
159 invD=1./D;
160 dvdX00=dXdv00*invD;
161 dvdX01=dXdv10*invD;
162 for (s=0;s<numTest; s++) {
163 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00;
164 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01;
165 }
166 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
167 }
168 }
169 }
170
171 }
172 #undef DIM
173 #undef LOCALDIM
174 }
175 /**************************************************************/
176 /* */
177 /* Jacobean 1D manifold in 2D and 2D elements */
178 /* */
179 void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
180 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
181 double* DSDv, dim_t numTest,double* DTDv,
182 double* dTdX, double* volume, index_t* element_id) {
183 #define DIM 2
184 #define LOCALDIM 2
185 register int e,q,s;
186 char error_msg[LenErrorMsg_MAX];
187 #pragma omp parallel
188 {
189 register double dXdv00,dXdv10,dXdv01,dXdv11,
190 dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
191 double X0[numShape], X1[numShape];
192 #pragma omp for private(e,q,s) schedule(static)
193 for(e=0;e<numElements;e++){
194 for (s=0;s<numShape; s++) {
195 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
196 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
197 }
198 for (q=0;q<numQuad;q++) {
199 dXdv00=0;
200 dXdv10=0;
201 dXdv01=0;
202 dXdv11=0;
203 for (s=0;s<numShape; s++) {
204 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
205 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
206 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
207 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
208 }
209 D = dXdv00*dXdv11 - dXdv01*dXdv10;
210 if (D==0.) {
211 sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
212 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
213 } else {
214 invD=1./D;
215 dvdX00= dXdv11*invD;
216 dvdX10=-dXdv10*invD;
217 dvdX01=-dXdv01*invD;
218 dvdX11= dXdv00*invD;
219
220 for (s=0;s<numTest; s++) {
221 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10;
222 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11;
223 }
224 }
225 volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
226 }
227 }
228
229 }
230 #undef DIM
231 #undef LOCALDIM
232 }
233 /**************************************************************/
234 /* */
235 /* Jacobean 3D */
236 /* */
237 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
238 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
239 double* DSDv, dim_t numTest,double* DTDv,
240 double* dTdX, double* volume, index_t* element_id) {
241 #define DIM 3
242 #define LOCALDIM 3
243 register int e,q,s;
244 char error_msg[LenErrorMsg_MAX];
245 #pragma omp parallel
246 {
247 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
248 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
249 double X0[numShape], X1[numShape], X2[numShape];
250 #pragma omp for private(e,q,s) schedule(static)
251 for(e=0;e<numElements;e++){
252 for (s=0;s<numShape; s++) {
253 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
254 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
255 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
256 }
257 for (q=0;q<numQuad;q++) {
258 dXdv00=0;
259 dXdv10=0;
260 dXdv20=0;
261 dXdv01=0;
262 dXdv11=0;
263 dXdv21=0;
264 dXdv02=0;
265 dXdv12=0;
266 dXdv22=0;
267 for (s=0;s<numShape; s++) {
268 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
269 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
270 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
271 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
272 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
273 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
274 dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
275 dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
276 dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
277 }
278 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
279 if (D==0.) {
280 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
281 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
282 } else {
283 invD=1./D;
284 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
285 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
286 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
287 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
288 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
289 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
290 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
291 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
292 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
293
294 for (s=0;s<numTest; s++) {
295 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX20;
296 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX21;
297 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX22;
298 }
299 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
300 }
301 }
302 }
303
304 }
305 #undef DIM
306 #undef LOCALDIM
307 }
308 /**************************************************************/
309 /* */
310 /* Jacobean 2D manifold in 3D with 3D elements */
311 /* */
312 void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
313 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
314 double* DSDv, dim_t numTest,double* DTDv,
315 double* dTdX, double* volume, index_t* element_id) {
316 #define DIM 3
317 #define LOCALDIM 3
318 register int e,q,s;
319 char error_msg[LenErrorMsg_MAX];
320 #pragma omp parallel
321 {
322 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
323 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
324 double X0[numShape], X1[numShape], X2[numShape], m0, m1, m2;
325 #pragma omp for private(e,q,s) schedule(static)
326 for(e=0;e<numElements;e++){
327 for (s=0;s<numShape; s++) {
328 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
329 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
330 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
331 }
332 for (q=0;q<numQuad;q++) {
333 dXdv00=0;
334 dXdv10=0;
335 dXdv20=0;
336 dXdv01=0;
337 dXdv11=0;
338 dXdv21=0;
339 dXdv02=0;
340 dXdv12=0;
341 dXdv22=0;
342 for (s=0;s<numShape; s++) {
343 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
344 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
345 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
346 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
347 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
348 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
349 dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
350 dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
351 dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
352 }
353 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
354 if (D==0.) {
355 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
356 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
357 } else {
358 invD=1./D;
359 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
360 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
361 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
362 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
363 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
364 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
365 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
366 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
367 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
368
369 for (s=0;s<numTest; s++) {
370 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX20;
371 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX21;
372 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCALDIM)]*dvdX22;
373 }
374 }
375 m0=dXdv10*dXdv21-dXdv20*dXdv11;
376 m1=dXdv20*dXdv01-dXdv00*dXdv21;
377 m2=dXdv00*dXdv11-dXdv10*dXdv01;
378 volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
379 }
380 }
381
382 }
383 #undef DIM
384 #undef LOCALDIM
385 }
386 /**************************************************************/
387 /* */
388 /* Jacobean 2D manifold in 3D with 2D elements */
389 /* */
390 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
391 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
392 double* DSDv, dim_t numTest,double* DTDv,
393 double* dTdX, double* volume, index_t* element_id) {
394 #define DIM 3
395 #define LOCALDIM 2
396 register int e,q,s;
397 char error_msg[LenErrorMsg_MAX];
398 #pragma omp parallel
399 {
400 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
401 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;
402 double X0[numShape], X1[numShape], X2[numShape];
403 #pragma omp for private(e,q,s) schedule(static)
404 for(e=0;e<numElements;e++){
405 for (s=0;s<numShape; s++) {
406 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
407 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
408 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
409 }
410 for (q=0;q<numQuad;q++) {
411 dXdv00=0;
412 dXdv10=0;
413 dXdv20=0;
414 dXdv01=0;
415 dXdv11=0;
416 dXdv21=0;
417 for (s=0;s<numShape; s++) {
418 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
419 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
420 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
421 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
422 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
423 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
424 }
425 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
426 m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
427 m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
428 D=m00*m11-m01*m01;
429 if (D==0.) {
430 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
431 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
432 } else {
433 invD=1./D;
434 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
435 dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
436 dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
437 dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
438 dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
439 dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
440 for (s=0;s<numTest; s++) {
441 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX10;
442 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX11;
443 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCALDIM)]*dvdX12;
444 }
445 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
446 }
447 }
448 }
449
450 }
451 #undef DIM
452 #undef LOCALDIM
453 }
454 /*
455 * $Log$
456 *
457 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26