/[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 3169 - (hide annotations)
Thu Sep 9 03:21:35 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 45750 byte(s)
Some moving

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26