/[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 853 - (hide annotations)
Wed Sep 20 05:56:36 2006 UTC (13 years, 1 month ago) by gross
File MIME type: text/plain
File size: 44566 byte(s)
some performance problems wit openmp fixed.
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 781 double* DSDv, dim_t numTest, double* DTDv,
36 gross 776 double* dTdX, double* volume, index_t* element_id) {
37     #define DIM 1
38 gross 781 #define LOCDIM 1
39 gross 776 register int e,q,s;
40     char error_msg[LenErrorMsg_MAX];
41     #pragma omp parallel
42     {
43 gross 853 register double D,invD, X0_loc;
44     #pragma omp for private(e,q,s,D,invD,X0_loc) schedule(static)
45 gross 776 for(e=0;e<numElements;e++){
46     for (q=0;q<numQuad;q++) {
47     D=0;
48 gross 853 for (s=0;s<numShape; s++) {
49     X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
50     D+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
51     }
52 gross 776 if (D==0.) {
53     sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);
54     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
55     } else {
56     invD=1./D;
57 gross 781 for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*invD;
58 gross 776 }
59 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
60 gross 776 }
61     }
62 jgs 82
63 gross 776 }
64     #undef DIM
65 gross 781 #undef LOCDIM
66 jgs 82
67 gross 776 }
68     /**************************************************************/
69     /* */
70 gross 777 /* Jacobean 2D with area element */
71 gross 776 /* */
72     void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
73 gross 777 dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
74 gross 781 double* DSDv, dim_t numTest, double* DTDv,
75 gross 776 double* dTdX, double* volume, index_t* element_id) {
76     #define DIM 2
77 gross 781 #define LOCDIM 2
78 gross 776 register int e,q,s;
79     char error_msg[LenErrorMsg_MAX];
80     #pragma omp parallel
81     {
82     register double dXdv00,dXdv10,dXdv01,dXdv11,
83 gross 853 dvdX00,dvdX10,dvdX01,dvdX11, D,invD,
84     X0_loc, X1_loc;
85     #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
86 gross 776 for(e=0;e<numElements;e++){
87     for (q=0;q<numQuad;q++) {
88     dXdv00=0;
89     dXdv10=0;
90     dXdv01=0;
91     dXdv11=0;
92     for (s=0;s<numShape; s++) {
93 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
94     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
95     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
96     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
97     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
98     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
99 gross 776 }
100     D = dXdv00*dXdv11 - dXdv01*dXdv10;
101     if (D==0.) {
102     sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
103     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
104     } else {
105     invD=1./D;
106     dvdX00= dXdv11*invD;
107     dvdX10=-dXdv10*invD;
108     dvdX01=-dXdv01*invD;
109     dvdX11= dXdv00*invD;
110 jgs 82
111 gross 776 for (s=0;s<numTest; s++) {
112 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;
113     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;
114 gross 776 }
115     }
116 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
117 gross 776 }
118     }
119 jgs 82
120 gross 776 }
121     #undef DIM
122 gross 781 #undef LOCDIM
123 gross 776 }
124     /**************************************************************/
125     /* */
126 gross 777 /* Jacobean 1D manifold in 2D and 1D elements */
127 gross 776 /* */
128 gross 777 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,double* QuadWeights,
129     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
130 gross 781 double* DSDv, dim_t numTest, double* DTDv,
131 gross 777 double* dTdX, double* volume, index_t* element_id) {
132     #define DIM 2
133 gross 781 #define LOCDIM 1
134 gross 776 register int e,q,s;
135     char error_msg[LenErrorMsg_MAX];
136     #pragma omp parallel
137     {
138 gross 853 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD,
139     X0_loc, X1_loc;
140     #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
141 gross 776 for(e=0;e<numElements;e++){
142     for (q=0;q<numQuad;q++) {
143     dXdv00=0;
144     dXdv10=0;
145 gross 777 for (s=0;s<numShape; s++) {
146 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
147     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
148     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
149     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
150 gross 777 }
151     D=dXdv00*dXdv00+dXdv10*dXdv10;
152     if (D==0.) {
153     sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
154     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
155     } else {
156     invD=1./D;
157     dvdX00=dXdv00*invD;
158     dvdX01=dXdv10*invD;
159     for (s=0;s<numTest; s++) {
160 gross 781 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00;
161     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01;
162 gross 777 }
163     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
164     }
165     }
166     }
167    
168     }
169     #undef DIM
170 gross 781 #undef LOCDIM
171 gross 777 }
172     /**************************************************************/
173     /* */
174 gross 781 /* Jacobean 1D manifold in 2D and 1D elements woth contact */
175     /* */
176     void Assemble_jacobeans_2D_M1D_E1D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
177     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
178     double* DSDv, dim_t numTest, double* DTDv,
179     double* dTdX, double* volume, index_t* element_id) {
180     #define DIM 2
181     #define LOCDIM 1
182     register int e,q,s;
183     char error_msg[LenErrorMsg_MAX];
184     #pragma omp parallel
185     {
186 gross 853 register double dXdv00_0,dXdv10_0,dvdX00_0,dvdX01_0,D_0,invD_0,
187     dXdv00_1,dXdv10_1,dvdX00_1,dvdX01_1,D_1,invD_1,
188     X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1;
189     #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dvdX00_0,dvdX01_0,D_0,invD_0,dXdv00_1,dXdv10_1,dvdX00_1,dvdX01_1,D_1,invD_1,X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1) schedule(static)
190 gross 781 for(e=0;e<numElements;e++){
191     for (q=0;q<numQuad;q++) {
192     dXdv00_0=0;
193     dXdv10_0=0;
194     dXdv00_1=0;
195     dXdv10_1=0;
196     for (s=0;s<numShape; s++) {
197 gross 853 X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s ,e,numNodes)],DIM)];
198     X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s ,e,numNodes)],DIM)];
199     X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
200     X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
201     dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
202     dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
203     dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
204     dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
205 gross 781 }
206     D_0=dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0;
207     D_1=dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1;
208     if (D_0 == 0.i || D_1 == 0.) {
209     sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
210     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
211     } else {
212     invD_0=1./D_0;
213     dvdX00_0=dXdv00_0*invD_0;
214     dvdX01_0=dXdv10_0*invD_0;
215     dvdX00_1=dXdv00_1*invD_1;
216     dvdX01_1=dXdv10_1*invD_1;
217     for (s=0;s<numTest; s++) {
218     dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0;
219     dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0;
220     dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1;
221     dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1;
222     }
223     volume[INDEX2(q,e,numQuad)]=(sqrt(D_0)+sqrt(D_1))/2.*QuadWeights[q];
224     }
225     }
226     }
227    
228     }
229     #undef DIM
230     #undef LOCDIM
231     }
232     /**************************************************************/
233     /* */
234 gross 777 /* Jacobean 1D manifold in 2D and 2D elements */
235     /* */
236     void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
237     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
238     double* DSDv, dim_t numTest,double* DTDv,
239     double* dTdX, double* volume, index_t* element_id) {
240     #define DIM 2
241 gross 781 #define LOCDIM 2
242 gross 777 register int e,q,s;
243     char error_msg[LenErrorMsg_MAX];
244     #pragma omp parallel
245     {
246     register double dXdv00,dXdv10,dXdv01,dXdv11,
247 gross 853 dvdX00,dvdX10,dvdX01,dvdX11, D,invD,
248     X0_loc, X1_loc;
249     #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
250 gross 777 for(e=0;e<numElements;e++){
251     for (q=0;q<numQuad;q++) {
252     dXdv00=0;
253     dXdv10=0;
254 gross 776 dXdv01=0;
255     dXdv11=0;
256     for (s=0;s<numShape; s++) {
257 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
258     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
259     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
260     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
261     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
262     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
263 gross 776 }
264 gross 777 D = dXdv00*dXdv11 - dXdv01*dXdv10;
265 gross 776 if (D==0.) {
266 gross 777 sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
267 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
268     } else {
269     invD=1./D;
270 gross 777 dvdX00= dXdv11*invD;
271     dvdX10=-dXdv10*invD;
272     dvdX01=-dXdv01*invD;
273     dvdX11= dXdv00*invD;
274 jgs 82
275 gross 776 for (s=0;s<numTest; s++) {
276 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;
277     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;
278 gross 776 }
279     }
280 gross 777 volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
281 gross 776 }
282     }
283 jgs 82
284 gross 776 }
285     #undef DIM
286 gross 781 #undef LOCDIM
287 gross 776 }
288     /**************************************************************/
289     /* */
290 gross 781 /* Jacobean 1D manifold in 2D and 2D elements with contact */
291     /* */
292     void Assemble_jacobeans_2D_M1D_E2D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
293     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
294     double* DSDv, dim_t numTest,double* DTDv,
295     double* dTdX, double* volume, index_t* element_id) {
296     #define DIM 2
297     #define LOCDIM 2
298     register int e,q,s;
299     char error_msg[LenErrorMsg_MAX];
300     #pragma omp parallel
301     {
302     register double dXdv00_0,dXdv10_0,dXdv01_0,dXdv11_0,dvdX00_0,dvdX10_0,dvdX01_0,dvdX11_0, D_0,invD_0,
303 gross 853 dXdv00_1,dXdv10_1,dXdv01_1,dXdv11_1,dvdX00_1,dvdX10_1,dvdX01_1,dvdX11_1, D_1,invD_1,
304     X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1;
305     #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dXdv01_0,dXdv11_0,dvdX00_0,dvdX10_0,dvdX01_0,dvdX11_0, D_0,invD_0,dXdv00_1,dXdv10_1,dXdv01_1,dXdv11_1,dvdX00_1,dvdX10_1,dvdX01_1,dvdX11_1, D_1,invD_1,X0_loc_0, X1_loc_0, X0_loc_1, X1_loc_1) schedule(static)
306 gross 781 for(e=0;e<numElements;e++){
307     for (q=0;q<numQuad;q++) {
308     dXdv00_0=0;
309     dXdv10_0=0;
310     dXdv01_0=0;
311     dXdv11_0=0;
312     dXdv00_1=0;
313     dXdv10_1=0;
314     dXdv01_1=0;
315     dXdv11_1=0;
316     for (s=0;s<numShape; s++) {
317 gross 853 X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s ,e,numNodes)],DIM)];
318     X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s ,e,numNodes)],DIM)];
319     X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
320     X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
321     dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
322     dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
323     dXdv01_0+=X0_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
324     dXdv11_0+=X1_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
325     dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
326     dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
327     dXdv01_1+=X0_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
328     dXdv11_1+=X1_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
329 gross 781 }
330     D_0 = dXdv00_0*dXdv11_0 - dXdv01_0*dXdv10_0;
331     D_1 = dXdv00_1*dXdv11_1 - dXdv01_1*dXdv10_1;
332     if ( (D_0 ==0.) || (D_1 ==0.) ) {
333     sprintf(error_msg,"Assemble_jacobeans_2D_E2D_C: element %d (id %d) has area zero.",e,element_id[e]);
334     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
335     } else {
336     invD_0=1./D_0;
337     dvdX00_0= dXdv11_0*invD_0;
338     dvdX10_0=-dXdv10_0*invD_0;
339     dvdX01_0=-dXdv01_0*invD_0;
340     dvdX11_0= dXdv00_0*invD_0;
341     invD_1=1./D_1;
342     dvdX00_1= dXdv11_1*invD_1;
343     dvdX10_1=-dXdv10_1*invD_1;
344     dvdX01_1=-dXdv01_1*invD_1;
345     dvdX11_1= dXdv00_1*invD_1;
346     for (s=0;s<numTest; s++) {
347     dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=
348     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0;
349     dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=
350     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0;
351     dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=
352     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1;
353     dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=
354     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1;
355     }
356     }
357     volume[INDEX2(q,e,numQuad)]=(sqrt(dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0)+sqrt(dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1))/2.*QuadWeights[q];
358     }
359     }
360    
361     }
362     #undef DIM
363     #undef LOCDIM
364     }
365     /**************************************************************/
366     /* */
367 gross 777 /* Jacobean 3D */
368 gross 776 /* */
369 gross 777 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
370     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
371 gross 781 double* DSDv, dim_t numTest, double* DTDv,
372 gross 777 double* dTdX, double* volume, index_t* element_id) {
373 gross 776 #define DIM 3
374 gross 781 #define LOCDIM 3
375 gross 853 int e,q,s;
376 gross 776 char error_msg[LenErrorMsg_MAX];
377     #pragma omp parallel
378     {
379 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
380 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
381     X0_loc,X1_loc,X2_loc;
382    
383     #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)
384 gross 776 for(e=0;e<numElements;e++){
385     for (q=0;q<numQuad;q++) {
386     dXdv00=0;
387     dXdv10=0;
388     dXdv20=0;
389     dXdv01=0;
390     dXdv11=0;
391     dXdv21=0;
392 gross 777 dXdv02=0;
393     dXdv12=0;
394     dXdv22=0;
395 gross 776 for (s=0;s<numShape; s++) {
396 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
397     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
398     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
399     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
400     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
401     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
402     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
403     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
404     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
405     dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
406     dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
407     dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
408 gross 776 }
409 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
410 gross 776 if (D==0.) {
411 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
412 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
413     } else {
414     invD=1./D;
415 gross 777 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
416     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
417     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
418     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
419     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
420     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
421     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
422     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
423     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
424    
425 gross 776 for (s=0;s<numTest; s++) {
426 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+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
427     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;
428     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;
429 gross 776 }
430 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
431 gross 776 }
432     }
433     }
434 jgs 82
435 gross 776 }
436     #undef DIM
437 gross 781 #undef LOCDIM
438 gross 776 }
439     /**************************************************************/
440     /* */
441 gross 777 /* Jacobean 2D manifold in 3D with 3D elements */
442 gross 776 /* */
443 gross 777 void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
444     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
445     double* DSDv, dim_t numTest,double* DTDv,
446     double* dTdX, double* volume, index_t* element_id) {
447 gross 776 #define DIM 3
448 gross 781 #define LOCDIM 3
449 gross 776 register int e,q,s;
450     char error_msg[LenErrorMsg_MAX];
451     #pragma omp parallel
452     {
453 gross 781 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,
454 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
455     X0_loc, X1_loc, X2_loc;
456     #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)
457 gross 776 for(e=0;e<numElements;e++){
458     for (q=0;q<numQuad;q++) {
459     dXdv00=0;
460     dXdv10=0;
461     dXdv20=0;
462 gross 777 dXdv01=0;
463     dXdv11=0;
464     dXdv21=0;
465     dXdv02=0;
466     dXdv12=0;
467     dXdv22=0;
468 gross 776 for (s=0;s<numShape; s++) {
469 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
470     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
471     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
472     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
473     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
474     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
475     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
476     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
477     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
478     dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
479     dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
480     dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
481 gross 776 }
482 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
483 gross 776 if (D==0.) {
484 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
485 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
486     } else {
487     invD=1./D;
488 gross 777 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
489     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
490     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
491     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
492     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
493     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
494     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
495     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
496     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
497    
498 gross 776 for (s=0;s<numTest; s++) {
499 gross 781 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=
500     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
501     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=
502     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
503     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=
504     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
505 gross 776 }
506     }
507 gross 777 m0=dXdv10*dXdv21-dXdv20*dXdv11;
508     m1=dXdv20*dXdv01-dXdv00*dXdv21;
509     m2=dXdv00*dXdv11-dXdv10*dXdv01;
510     volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
511 gross 776 }
512 jgs 82 }
513    
514 gross 776 }
515     #undef DIM
516 gross 781 #undef LOCDIM
517 gross 776 }
518     /**************************************************************/
519     /* */
520 gross 781 /* Jacobean 2D manifold in 3D with 3D elements on contact */
521     /* */
522     void Assemble_jacobeans_3D_M2D_E3D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
523     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
524     double* DSDv, dim_t numTest,double* DTDv,
525     double* dTdX, double* volume, index_t* element_id) {
526     #define DIM 3
527     #define LOCDIM 3
528     register int e,q,s;
529     char error_msg[LenErrorMsg_MAX];
530     #pragma omp parallel
531     {
532     register double dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,dXdv02_0,dXdv12_0,dXdv22_0, m0_0, m1_0, m2_0,
533     dvdX00_0,dvdX10_0,dvdX20_0,dvdX01_0,dvdX11_0,dvdX21_0,dvdX02_0,dvdX12_0,dvdX22_0, D_0,invD_0,
534     dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,dXdv02_1,dXdv12_1,dXdv22_1, m0_1, m1_1, m2_1,
535 gross 853 dvdX00_1,dvdX10_1,dvdX20_1,dvdX01_1,dvdX11_1,dvdX21_1,dvdX02_1,dvdX12_1,dvdX22_1, D_1,invD_1,
536     X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1;
537     #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,dXdv02_0,dXdv12_0,dXdv22_0, m0_0, m1_0, m2_0,dvdX00_0,dvdX10_0,dvdX20_0,dvdX01_0,dvdX11_0,dvdX21_0,dvdX02_0,dvdX12_0,dvdX22_0, D_0,invD_0,dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,dXdv02_1,dXdv12_1,dXdv22_1, m0_1, m1_1, m2_1,dvdX00_1,dvdX10_1,dvdX20_1,dvdX01_1,dvdX11_1,dvdX21_1,dvdX02_1,dvdX12_1,dvdX22_1, D_1,invD_1,X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1) schedule(static)
538 gross 781 for(e=0;e<numElements;e++){
539     for (q=0;q<numQuad;q++) {
540     dXdv00_0=0;
541     dXdv10_0=0;
542     dXdv20_0=0;
543     dXdv01_0=0;
544     dXdv11_0=0;
545     dXdv21_0=0;
546     dXdv02_0=0;
547     dXdv12_0=0;
548     dXdv22_0=0;
549     dXdv00_1=0;
550     dXdv10_1=0;
551     dXdv20_1=0;
552     dXdv01_1=0;
553     dXdv11_1=0;
554     dXdv21_1=0;
555     dXdv02_1=0;
556     dXdv12_1=0;
557     dXdv22_1=0;
558     for (s=0;s<numShape; s++) {
559 gross 853 X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s ,e,numNodes)],DIM)];
560     X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s ,e,numNodes)],DIM)];
561     X2_loc_0=coordinates[INDEX2(2,nodes[INDEX2(s ,e,numNodes)],DIM)];
562     X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
563     X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
564     X2_loc_1=coordinates[INDEX2(2,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
565     dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
566     dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
567     dXdv20_0+=X2_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
568     dXdv01_0+=X0_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
569     dXdv11_0+=X1_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
570     dXdv21_0+=X2_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
571     dXdv02_0+=X0_loc_0*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
572     dXdv12_0+=X1_loc_0*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
573     dXdv22_0+=X2_loc_0*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
574     dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
575     dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
576     dXdv20_1+=X2_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
577     dXdv01_1+=X0_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
578     dXdv11_1+=X1_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
579     dXdv21_1+=X2_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
580     dXdv02_1+=X0_loc_1*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
581     dXdv12_1+=X1_loc_1*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
582     dXdv22_1+=X2_loc_1*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
583 gross 781 }
584    
585     D_0=dXdv00_0*(dXdv11_0*dXdv22_0-dXdv12_0*dXdv21_0)+dXdv01_0*(dXdv20_0*dXdv12_0-dXdv10_0*dXdv22_0)+dXdv02_0*(dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0);
586     D_1=dXdv00_1*(dXdv11_1*dXdv22_1-dXdv12_1*dXdv21_1)+dXdv01_1*(dXdv20_1*dXdv12_1-dXdv10_1*dXdv22_1)+dXdv02_1*(dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1);
587     if ( (D_0==0.) || (D_1 == 0.)) {
588     sprintf(error_msg,"Assemble_jacobeans_3D_C: element %d (id %d) has volume zero.",e,element_id[e]);
589     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
590     } else {
591     invD_0=1./D_0;
592     dvdX00_0=(dXdv11_0*dXdv22_0-dXdv12_0*dXdv21_0)*invD_0;
593     dvdX10_0=(dXdv20_0*dXdv12_0-dXdv10_0*dXdv22_0)*invD_0;
594     dvdX20_0=(dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0)*invD_0;
595     dvdX01_0=(dXdv02_0*dXdv21_0-dXdv01_0*dXdv22_0)*invD_0;
596     dvdX11_0=(dXdv00_0*dXdv22_0-dXdv20_0*dXdv02_0)*invD_0;
597     dvdX21_0=(dXdv01_0*dXdv20_0-dXdv00_0*dXdv21_0)*invD_0;
598     dvdX02_0=(dXdv01_0*dXdv12_0-dXdv02_0*dXdv11_0)*invD_0;
599     dvdX12_0=(dXdv02_0*dXdv10_0-dXdv00_0*dXdv12_0)*invD_0;
600     dvdX22_0=(dXdv00_0*dXdv11_0-dXdv01_0*dXdv10_0)*invD_0;
601     invD_1=1./D_1;
602     dvdX00_1=(dXdv11_1*dXdv22_1-dXdv12_1*dXdv21_1)*invD_1;
603     dvdX10_1=(dXdv20_1*dXdv12_1-dXdv10_1*dXdv22_1)*invD_1;
604     dvdX20_1=(dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1)*invD_1;
605     dvdX01_1=(dXdv02_1*dXdv21_1-dXdv01_1*dXdv22_1)*invD_1;
606     dvdX11_1=(dXdv00_1*dXdv22_1-dXdv20_1*dXdv02_1)*invD_1;
607     dvdX21_1=(dXdv01_1*dXdv20_1-dXdv00_1*dXdv21_1)*invD_1;
608     dvdX02_1=(dXdv01_1*dXdv12_1-dXdv02_1*dXdv11_1)*invD_1;
609     dvdX12_1=(dXdv02_1*dXdv10_1-dXdv00_1*dXdv12_1)*invD_1;
610     dvdX22_1=(dXdv00_1*dXdv11_1-dXdv01_1*dXdv10_1)*invD_1;
611    
612     for (s=0;s<numTest; s++) {
613     dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=
614     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20_0;
615     dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=
616     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21_0;
617     dTdX[INDEX4( s,2,q,e,2*numTest,DIM,numQuad)]=
618     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_0+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22_0;
619     dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=
620     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20_1;
621     dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=
622     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21_1;
623     dTdX[INDEX4(numTest+s,2,q,e,2*numTest,DIM,numQuad)]=
624     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_1+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22_1;
625     }
626     }
627     m0_0=dXdv10_0*dXdv21_0-dXdv20_0*dXdv11_0;
628     m1_0=dXdv20_0*dXdv01_0-dXdv00_0*dXdv21_0;
629     m2_0=dXdv00_0*dXdv11_0-dXdv10_0*dXdv01_0;
630     m0_1=dXdv10_1*dXdv21_1-dXdv20_1*dXdv11_1;
631     m1_1=dXdv20_1*dXdv01_1-dXdv00_1*dXdv21_1;
632     m2_1=dXdv00_1*dXdv11_1-dXdv10_1*dXdv01_1;
633     volume[INDEX2(q,e,numQuad)]=(sqrt(m0_0*m0_0+m1_0*m1_0+m2_0*m2_0)+sqrt(m0_1*m0_1+m1_1*m1_1+m2_1*m2_1))/2.*QuadWeights[q];
634     }
635     }
636    
637     }
638     #undef DIM
639     #undef LOCDIM
640     }
641     /**************************************************************/
642     /* */
643 gross 777 /* Jacobean 2D manifold in 3D with 2D elements */
644 gross 776 /* */
645 gross 777 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
646     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
647     double* DSDv, dim_t numTest,double* DTDv,
648     double* dTdX, double* volume, index_t* element_id) {
649     #define DIM 3
650 gross 781 #define LOCDIM 2
651 gross 776 register int e,q,s;
652     char error_msg[LenErrorMsg_MAX];
653     #pragma omp parallel
654 jgs 82 {
655 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
656 gross 853 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
657     X0_loc, X1_loc, X2_loc;
658     #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)
659 gross 776 for(e=0;e<numElements;e++){
660     for (q=0;q<numQuad;q++) {
661     dXdv00=0;
662     dXdv10=0;
663 gross 777 dXdv20=0;
664     dXdv01=0;
665     dXdv11=0;
666     dXdv21=0;
667 gross 776 for (s=0;s<numShape; s++) {
668 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
669     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
670     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
671     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
672     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
673     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
674     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
675     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
676     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
677 gross 776 }
678 gross 777 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
679     m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
680     m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
681     D=m00*m11-m01*m01;
682 gross 776 if (D==0.) {
683 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
684 gross 776 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
685     } else {
686     invD=1./D;
687 gross 777 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
688     dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
689     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
690     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
691     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
692     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
693 gross 776 for (s=0;s<numTest; s++) {
694 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;
695     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;
696     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;
697 gross 776 }
698     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
699     }
700     }
701     }
702 jgs 82
703     }
704 gross 776 #undef DIM
705 gross 781 #undef LOCDIM
706 jgs 82 }
707 gross 781 /**************************************************************/
708     /* */
709     /* Jacobean 2D manifold in 3D with 2D elements with contact */
710     /* */
711     void Assemble_jacobeans_3D_M2D_E2D_C(double* coordinates, dim_t numQuad,double* QuadWeights,
712     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
713     double* DSDv, dim_t numTest,double* DTDv,
714     double* dTdX, double* volume, index_t* element_id) {
715     #define DIM 3
716     #define LOCDIM 2
717     register int e,q,s;
718     char error_msg[LenErrorMsg_MAX];
719     #pragma omp parallel
720     {
721     register double dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,m00_0,m01_0,m11_0,
722     dvdX00_0,dvdX01_0,dvdX02_0,dvdX10_0,dvdX11_0,dvdX12_0,D_0,invD_0,
723     dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,m00_1,m01_1,m11_1,
724 gross 853 dvdX00_1,dvdX01_1,dvdX02_1,dvdX10_1,dvdX11_1,dvdX12_1,D_1,invD_1,
725     X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1;
726     #pragma omp for private(e,q,s,dXdv00_0,dXdv10_0,dXdv20_0,dXdv01_0,dXdv11_0,dXdv21_0,m00_0,m01_0,m11_0,dvdX00_0,dvdX01_0,dvdX02_0,dvdX10_0,dvdX11_0,dvdX12_0,D_0,invD_0,dXdv00_1,dXdv10_1,dXdv20_1,dXdv01_1,dXdv11_1,dXdv21_1,m00_1,m01_1,m11_1,dvdX00_1,dvdX01_1,dvdX02_1,dvdX10_1,dvdX11_1,dvdX12_1,D_1,invD_1,X0_loc_0, X1_loc_0, X2_loc_0, X0_loc_1, X1_loc_1, X2_loc_1) schedule(static)
727 gross 781 for(e=0;e<numElements;e++){
728     for (q=0;q<numQuad;q++) {
729     dXdv00_0=0;
730     dXdv10_0=0;
731     dXdv20_0=0;
732     dXdv01_0=0;
733     dXdv11_0=0;
734     dXdv21_0=0;
735     dXdv00_1=0;
736     dXdv10_1=0;
737     dXdv20_1=0;
738     dXdv01_1=0;
739     dXdv11_1=0;
740     dXdv21_1=0;
741     for (s=0;s<numShape; s++) {
742 gross 853 X0_loc_0=coordinates[INDEX2(0,nodes[INDEX2(s ,e,numNodes)],DIM)];
743     X1_loc_0=coordinates[INDEX2(1,nodes[INDEX2(s ,e,numNodes)],DIM)];
744     X2_loc_0=coordinates[INDEX2(2,nodes[INDEX2(s ,e,numNodes)],DIM)];
745     X0_loc_1=coordinates[INDEX2(0,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
746     X1_loc_1=coordinates[INDEX2(1,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
747     X2_loc_1=coordinates[INDEX2(2,nodes[INDEX2(s+numShape,e,numNodes)],DIM)];
748     dXdv00_0+=X0_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
749     dXdv10_0+=X1_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
750     dXdv20_0+=X2_loc_0*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
751     dXdv01_0+=X0_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
752     dXdv11_0+=X1_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
753     dXdv21_0+=X2_loc_0*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
754     dXdv00_1+=X0_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
755     dXdv10_1+=X1_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
756     dXdv20_1+=X2_loc_1*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
757     dXdv01_1+=X0_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
758     dXdv11_1+=X1_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
759     dXdv21_1+=X2_loc_1*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
760 gross 781 }
761     m00_0=dXdv00_0*dXdv00_0+dXdv10_0*dXdv10_0+dXdv20_0*dXdv20_0;
762     m01_0=dXdv00_0*dXdv01_0+dXdv10_0*dXdv11_0+dXdv20_0*dXdv21_0;
763     m11_0=dXdv01_0*dXdv01_0+dXdv11_0*dXdv11_0+dXdv21_0*dXdv21_0;
764     D_0=m00_0*m11_0-m01_0*m01_0;
765     m00_1=dXdv00_1*dXdv00_1+dXdv10_1*dXdv10_1+dXdv20_1*dXdv20_1;
766     m01_1=dXdv00_1*dXdv01_1+dXdv10_1*dXdv11_1+dXdv20_1*dXdv21_1;
767     m11_1=dXdv01_1*dXdv01_1+dXdv11_1*dXdv11_1+dXdv21_1*dXdv21_1;
768     D_1=m00_1*m11_1-m01_1*m01_1;
769     if ( (D_0==0.) || (D_1 == 0.) ) {
770     sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
771     Finley_setError(ZERO_DIVISION_ERROR,error_msg);
772     } else {
773     invD_0=1./D_0;
774     dvdX00_0=( m00_0*dXdv00_0-m01_0*dXdv01_0)*invD_0;
775     dvdX01_0=( m00_0*dXdv10_0-m01_0*dXdv11_0)*invD_0;
776     dvdX02_0=( m00_0*dXdv20_0-m01_0*dXdv21_0)*invD_0;
777     dvdX10_0=(-m01_0*dXdv00_0+m11_0*dXdv01_0)*invD_0;
778     dvdX11_0=(-m01_0*dXdv10_0+m11_0*dXdv11_0)*invD_0;
779     dvdX12_0=(-m01_0*dXdv20_0+m11_0*dXdv21_0)*invD_0;
780     invD_1=1./D_1;
781     dvdX00_1=( m00_1*dXdv00_1-m01_1*dXdv01_1)*invD_1;
782     dvdX01_1=( m00_1*dXdv10_1-m01_1*dXdv11_1)*invD_1;
783     dvdX02_1=( m00_1*dXdv20_1-m01_1*dXdv21_1)*invD_1;
784     dvdX10_1=(-m01_1*dXdv00_1+m11_1*dXdv01_1)*invD_1;
785     dvdX11_1=(-m01_1*dXdv10_1+m11_1*dXdv11_1)*invD_1;
786     dvdX12_1=(-m01_1*dXdv20_1+m11_1*dXdv21_1)*invD_1;
787     for (s=0;s<numTest; s++) {
788     dTdX[INDEX4( s,0,q,e,2*numTest,DIM,numQuad)]=
789     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_0;
790     dTdX[INDEX4( s,1,q,e,2*numTest,DIM,numQuad)]=
791     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_0;
792     dTdX[INDEX4( s,2,q,e,2*numTest,DIM,numQuad)]=
793     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_0+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_0;
794     dTdX[INDEX4(numTest+s,0,q,e,2*numTest,DIM,numQuad)]=
795     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10_1;
796     dTdX[INDEX4(numTest+s,1,q,e,2*numTest,DIM,numQuad)]=
797     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11_1;
798     dTdX[INDEX4(numTest+s,2,q,e,2*numTest,DIM,numQuad)]=
799     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02_1+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12_1;
800     }
801     volume[INDEX2(q,e,numQuad)]=(sqrt(D_0)+sqrt(D_1))/2.*QuadWeights[q];
802     }
803     }
804     }
805    
806     }
807     #undef DIM
808     #undef LOCDIM
809     }
810 jgs 82 /*
811 gross 853 * $Log:$
812 jgs 150 *
813 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26