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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26