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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26