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

Annotation of /trunk/finley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 777 - (hide annotations)
Wed Jul 12 08:54:45 2006 UTC (13 years, 3 months ago) by gross
File MIME type: text/plain
File size: 21708 byte(s)
integration with persistent jacobeans is running now
1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
15    
16 jgs 150 /* Author: gross@access.edu.au */
17     /* Version: $Id$ */
18 jgs 82
19     /**************************************************************/
20    
21     #include "Assemble.h"
22     #include "Util.h"
23     #ifdef _OPENMP
24     #include <omp.h>
25     #endif
26    
27 gross 777
28    
29 jgs 82 /**************************************************************/
30 gross 776 /* */
31     /* Jacobean 1D */
32     /* */
33     void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
34 gross 777 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
35 gross 776 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 gross 777 for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
48 gross 776 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 gross 777 for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;
57 gross 776 }
58 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
59 gross 776 }
60     }
61 jgs 82
62 gross 776 }
63     #undef DIM
64     #undef LOCALDIM
65 jgs 82
66 gross 776 }
67     /**************************************************************/
68     /* */
69 gross 777 /* Jacobean 2D with area element */
70 gross 776 /* */
71     void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
72 gross 777 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
73 gross 776 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 gross 777 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
88     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
89 gross 776 }
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 jgs 82
112 gross 776 for (s=0;s<numTest; s++) {
113 gross 777 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 gross 776 }
116     }
117 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
118 gross 776 }
119     }
120 jgs 82
121 gross 776 }
122     #undef DIM
123     #undef LOCALDIM
124     }
125     /**************************************************************/
126     /* */
127 gross 777 /* Jacobean 1D manifold in 2D and 1D elements */
128 gross 776 /* */
129 gross 777 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 gross 776 register int e,q,s;
136     char error_msg[LenErrorMsg_MAX];
137     #pragma omp parallel
138     {
139 gross 777 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
140     double X0[numShape], X1[numShape];
141 gross 776 #pragma omp for private(e,q,s) schedule(static)
142     for(e=0;e<numElements;e++){
143     for (s=0;s<numShape; s++) {
144 gross 777 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
145     X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
146 gross 776 }
147     for (q=0;q<numQuad;q++) {
148     dXdv00=0;
149     dXdv10=0;
150 gross 777 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 gross 776 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 gross 777 D = dXdv00*dXdv11 - dXdv01*dXdv10;
210 gross 776 if (D==0.) {
211 gross 777 sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
212 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
213     } else {
214     invD=1./D;
215 gross 777 dvdX00= dXdv11*invD;
216     dvdX10=-dXdv10*invD;
217     dvdX01=-dXdv01*invD;
218     dvdX11= dXdv00*invD;
219 jgs 82
220 gross 776 for (s=0;s<numTest; s++) {
221 gross 777 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 gross 776 }
224     }
225 gross 777 volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
226 gross 776 }
227     }
228 jgs 82
229 gross 776 }
230     #undef DIM
231     #undef LOCALDIM
232     }
233     /**************************************************************/
234     /* */
235 gross 777 /* Jacobean 3D */
236 gross 776 /* */
237 gross 777 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 gross 776 #define DIM 3
242 gross 777 #define LOCALDIM 3
243 gross 776 register int e,q,s;
244     char error_msg[LenErrorMsg_MAX];
245     #pragma omp parallel
246     {
247 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
248     dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
249 gross 776 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 gross 777 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 gross 776 }
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 gross 777 dXdv02=0;
265     dXdv12=0;
266     dXdv22=0;
267 gross 776 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 gross 777 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 gross 776 }
278 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
279 gross 776 if (D==0.) {
280 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
281 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
282     } else {
283     invD=1./D;
284 gross 777 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 gross 776 for (s=0;s<numTest; s++) {
295 gross 777 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 gross 776 }
299 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
300 gross 776 }
301     }
302     }
303 jgs 82
304 gross 776 }
305     #undef DIM
306     #undef LOCALDIM
307     }
308     /**************************************************************/
309     /* */
310 gross 777 /* Jacobean 2D manifold in 3D with 3D elements */
311 gross 776 /* */
312 gross 777 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 gross 776 #define DIM 3
317 gross 777 #define LOCALDIM 3
318 gross 776 register int e,q,s;
319     char error_msg[LenErrorMsg_MAX];
320     #pragma omp parallel
321     {
322 gross 777 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 gross 776 #pragma omp for private(e,q,s) schedule(static)
326     for(e=0;e<numElements;e++){
327     for (s=0;s<numShape; s++) {
328 gross 777 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 gross 776 }
332     for (q=0;q<numQuad;q++) {
333     dXdv00=0;
334     dXdv10=0;
335     dXdv20=0;
336 gross 777 dXdv01=0;
337     dXdv11=0;
338     dXdv21=0;
339     dXdv02=0;
340     dXdv12=0;
341     dXdv22=0;
342 gross 776 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 gross 777 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 gross 776 }
353 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
354 gross 776 if (D==0.) {
355 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
356 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
357     } else {
358     invD=1./D;
359 gross 777 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 gross 776 for (s=0;s<numTest; s++) {
370 gross 777 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 gross 776 }
374     }
375 gross 777 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 gross 776 }
380 jgs 82 }
381    
382 gross 776 }
383     #undef DIM
384     #undef LOCALDIM
385     }
386     /**************************************************************/
387     /* */
388 gross 777 /* Jacobean 2D manifold in 3D with 2D elements */
389 gross 776 /* */
390 gross 777 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 gross 776 register int e,q,s;
397     char error_msg[LenErrorMsg_MAX];
398     #pragma omp parallel
399 jgs 82 {
400 gross 777 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 gross 776 #pragma omp for private(e,q,s) schedule(static)
404     for(e=0;e<numElements;e++){
405     for (s=0;s<numShape; s++) {
406 gross 777 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 gross 776 }
410     for (q=0;q<numQuad;q++) {
411     dXdv00=0;
412     dXdv10=0;
413 gross 777 dXdv20=0;
414     dXdv01=0;
415     dXdv11=0;
416     dXdv21=0;
417 gross 776 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 gross 777 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 gross 776 }
425 gross 777 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 gross 776 if (D==0.) {
430 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
431 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
432     } else {
433     invD=1./D;
434 gross 777 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 gross 776 for (s=0;s<numTest; s++) {
441 gross 777 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 gross 776 }
445     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
446     }
447     }
448     }
449 jgs 82
450     }
451 gross 776 #undef DIM
452     #undef LOCALDIM
453 jgs 82 }
454     /*
455     * $Log$
456 jgs 150 *
457 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26