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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2750 - (hide annotations)
Tue Nov 17 23:39:20 2009 UTC (9 years, 6 months ago) by caltinay
Original Path: trunk/finley/src/Assemble_jacobeans.c
File MIME type: text/plain
File size: 44733 byte(s)
Fixed some compilation errors with OpenMP introduced in last commit.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26