/[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 3171 - (hide annotations)
Fri Sep 10 00:31:11 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 20970 byte(s)
Remove some contact jacobean functions.
Assemble_jacobeans_2D_M1D_E2D finished first pass

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 gross 777 /* Jacobean 2D with area element */
48 gross 776 /* */
49 jfenwick 3170 void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
50     double* dTdX, double* volume, index_t* element_id)
51     {
52     #define DIM 2
53     #define LOCDIM 2
54     register int e,q,s;
55     char error_msg[LenErrorMsg_MAX];
56     const double quadweight=(numQuad==1)?1./2:1./6; /* numQuad is 1 or 3 */
57     const dim_t numTest=3; // hoping this is used in constant folding
58     static const int DTDV[3][2]={{-1,-1}, {1,0}, {0,1}}; // if constant folding is not applied here will need to write
59     // out in full
60     #pragma omp parallel
61     {
62     register double dXdv00,dXdv10,dXdv01,dXdv11,
63     dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
64     #pragma omp for private(e,q,s,dXdv00,dXdv10,dXdv01,dXdv11,dvdX00,dvdX10,dvdX01,dvdX11, D,invD,X0_loc, X1_loc) schedule(static)
65     for(e=0;e<numElements;e++)
66     {
67 jfenwick 3169 #define COMPDXDV0(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
68     coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*1+\
69     coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(0)
70    
71     #define COMPDXDV1(P) coordinates[INDEX2(P,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1)+\
72     coordinates[INDEX2(P,nodes[INDEX2(1,e,numNodes)],DIM)]*(0)+\
73     coordinates[INDEX2(P,nodes[INDEX2(2,e,numNodes)],DIM)]*(1)
74    
75 jfenwick 3170 dXdv00=0;
76     dXdv10=0;
77     dXdv01=0;
78     dXdv11=0;
79     dXdv00=COMPDXDV0(0);
80     dXdv10=COMPDXDV0(1);
81     dXdv01=COMPDXDV1(0);
82     dXdv11=COMPDXDV1(1);
83     D = dXdv00*dXdv11 - dXdv01*dXdv10;
84     if (D==0.)
85     {
86     sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
87     Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
88     }
89     else
90     {
91     invD=1./D;
92     dvdX00= dXdv11*invD;
93     dvdX10=-dXdv10*invD;
94     dvdX01=-dXdv01*invD;
95     dvdX11= dXdv00*invD;
96     if (numQuad==1)
97     {
98     for (s=0;s<3; s++)
99 jfenwick 3169 {
100 jfenwick 3170 #define DTDXSET(P,Q) dTdX[INDEX4(s,P,Q,e,numTest,DIM,numQuad)] = DTDV[s][0]*dvdX0##P+DTDV[s][1]*dvdX1##P
101    
102     DTDXSET(0,0);
103     DTDXSET(1,0);
104     }
105     volume[INDEX2(0,e,1)]=ABS(D)*quadweight;
106     }
107     else // numQuad==3
108     {
109     for (q=0;q<numTest;++q) // relying on unroll loops to optimise this
110     {
111     for (s=0;s<3; s++)
112     {
113     DTDXSET(0,q);
114     DTDXSET(1,q);
115 jfenwick 3169 }
116 jfenwick 3170 volume[INDEX2(q,e,3)]=ABS(D)*quadweight;
117 jfenwick 3169 }
118     }
119 gross 776 }
120 jfenwick 3169 }
121 gross 776 }
122     #undef DIM
123 gross 781 #undef LOCDIM
124 jfenwick 3170 #undef DTDXSET
125 jfenwick 3169 #undef COMPDXDV0
126     #undef COMPDXDV1
127 gross 776 }
128     /**************************************************************/
129     /* */
130 gross 777 /* Jacobean 1D manifold in 2D and 1D elements */
131 gross 776 /* */
132 jfenwick 3171 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,/*double* bogus_QuadWeights,
133     dim_t bogus_numShape,*/ dim_t numElements, dim_t numNodes, index_t* nodes,
134     /*double* bogus_DSDv, dim_t bogus_numTest, double* DTDv,*/
135 gross 777 double* dTdX, double* volume, index_t* element_id) {
136     #define DIM 2
137 gross 781 #define LOCDIM 1
138 jfenwick 3171 register int e;
139 gross 776 char error_msg[LenErrorMsg_MAX];
140 jfenwick 3171 const dim_t numTest=2;
141     double quadweight=(numQuad==1)?1.0:0.5;
142     // numQuad is 1 or 2
143 gross 776 #pragma omp parallel
144     {
145 jfenwick 3171 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
146 caltinay 2750 #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
147 gross 776 for(e=0;e<numElements;e++){
148     dXdv00=0;
149     dXdv10=0;
150 jfenwick 3171 dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
151     dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
152 gross 777 D=dXdv00*dXdv00+dXdv10*dXdv10;
153 jfenwick 3171 if (D==0.) {
154 gross 777 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
155 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
156 jfenwick 3171 } else {
157 gross 777 invD=1./D;
158     dvdX00=dXdv00*invD;
159     dvdX01=dXdv10*invD;
160 jfenwick 3171 // The number of quad points is 1 or 2
161 gross 777
162 jfenwick 3170
163 jfenwick 3171 dTdX[INDEX4(0,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
164     dTdX[INDEX4(0,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
165     dTdX[INDEX4(1,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
166     dTdX[INDEX4(1,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
167     volume[INDEX2(0,e,numQuad)]=sqrt(D)*quadweight;
168    
169     if (numQuad==2)
170     {
171    
172     dTdX[INDEX4(0,0,1,e,numTest,DIM,numQuad)]=dvdX00;
173     dTdX[INDEX4(0,1,1,e,numTest,DIM,numQuad)]=dvdX01;
174     dTdX[INDEX4(1,0,1,e,numTest,DIM,numQuad)]=dvdX00;
175     dTdX[INDEX4(1,1,1,e,numTest,DIM,numQuad)]=dvdX01;
176     volume[INDEX2(1,e,numQuad)]=sqrt(D)*quadweight;
177    
178     }
179 gross 781 }
180     }
181    
182     }
183     #undef DIM
184     #undef LOCDIM
185     }
186 jfenwick 3171
187 gross 781 /**************************************************************/
188     /* */
189 gross 777 /* Jacobean 1D manifold in 2D and 2D elements */
190     /* */
191     void Assemble_jacobeans_2D_M1D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
192     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
193     double* DSDv, dim_t numTest,double* DTDv,
194     double* dTdX, double* volume, index_t* element_id) {
195     #define DIM 2
196 gross 781 #define LOCDIM 2
197 gross 777 register int e,q,s;
198     char error_msg[LenErrorMsg_MAX];
199     #pragma omp parallel
200     {
201     register double dXdv00,dXdv10,dXdv01,dXdv11,
202 gross 853 dvdX00,dvdX10,dvdX01,dvdX11, D,invD,
203     X0_loc, X1_loc;
204 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)
205 gross 777 for(e=0;e<numElements;e++){
206     for (q=0;q<numQuad;q++) {
207     dXdv00=0;
208     dXdv10=0;
209 gross 776 dXdv01=0;
210     dXdv11=0;
211     for (s=0;s<numShape; s++) {
212 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
213     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
214     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
215     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
216     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
217     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
218 gross 776 }
219 gross 777 D = dXdv00*dXdv11 - dXdv01*dXdv10;
220 gross 776 if (D==0.) {
221 gross 777 sprintf(error_msg,"Assemble_jacobeans_2D_E2D: element %d (id %d) has area zero.",e,element_id[e]);
222 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
223 gross 776 } else {
224     invD=1./D;
225 gross 777 dvdX00= dXdv11*invD;
226     dvdX10=-dXdv10*invD;
227     dvdX01=-dXdv01*invD;
228     dvdX11= dXdv00*invD;
229 jgs 82
230 gross 776 for (s=0;s<numTest; s++) {
231 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;
232     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;
233 gross 776 }
234     }
235 gross 777 volume[INDEX2(q,e,numQuad)]=sqrt(dXdv00*dXdv00+dXdv10*dXdv10)*QuadWeights[q];
236 gross 776 }
237     }
238 jgs 82
239 gross 776 }
240     #undef DIM
241 gross 781 #undef LOCDIM
242 gross 776 }
243 gross 781
244     /**************************************************************/
245     /* */
246 gross 777 /* Jacobean 3D */
247 gross 776 /* */
248 gross 777 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
249     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
250 gross 781 double* DSDv, dim_t numTest, double* DTDv,
251 gross 777 double* dTdX, double* volume, index_t* element_id) {
252 gross 776 #define DIM 3
253 gross 781 #define LOCDIM 3
254 gross 853 int e,q,s;
255 gross 776 char error_msg[LenErrorMsg_MAX];
256     #pragma omp parallel
257     {
258 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
259 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
260     X0_loc,X1_loc,X2_loc;
261    
262 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)
263 gross 776 for(e=0;e<numElements;e++){
264     for (q=0;q<numQuad;q++) {
265     dXdv00=0;
266     dXdv10=0;
267     dXdv20=0;
268     dXdv01=0;
269     dXdv11=0;
270     dXdv21=0;
271 gross 777 dXdv02=0;
272     dXdv12=0;
273     dXdv22=0;
274 gross 776 for (s=0;s<numShape; s++) {
275 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
276     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
277     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
278     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
279     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
280     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
281     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
282     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
283     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
284     dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
285     dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
286     dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
287 gross 776 }
288 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
289 gross 776 if (D==0.) {
290 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
291 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
292 gross 776 } else {
293     invD=1./D;
294 gross 777 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
295     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
296     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
297     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
298     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
299     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
300     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
301     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
302     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
303    
304 gross 776 for (s=0;s<numTest; s++) {
305 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;
306     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;
307     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;
308 gross 776 }
309 gross 777 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
310 gross 776 }
311     }
312     }
313 jgs 82
314 gross 776 }
315     #undef DIM
316 gross 781 #undef LOCDIM
317 gross 776 }
318     /**************************************************************/
319     /* */
320 gross 777 /* Jacobean 2D manifold in 3D with 3D elements */
321 gross 776 /* */
322 gross 777 void Assemble_jacobeans_3D_M2D_E3D(double* coordinates, dim_t numQuad,double* QuadWeights,
323     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
324     double* DSDv, dim_t numTest,double* DTDv,
325     double* dTdX, double* volume, index_t* element_id) {
326 gross 776 #define DIM 3
327 gross 781 #define LOCDIM 3
328 gross 776 register int e,q,s;
329     char error_msg[LenErrorMsg_MAX];
330     #pragma omp parallel
331     {
332 gross 781 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22, m0, m1, m2,
333 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
334     X0_loc, X1_loc, X2_loc;
335 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)
336 gross 776 for(e=0;e<numElements;e++){
337     for (q=0;q<numQuad;q++) {
338     dXdv00=0;
339     dXdv10=0;
340     dXdv20=0;
341 gross 777 dXdv01=0;
342     dXdv11=0;
343     dXdv21=0;
344     dXdv02=0;
345     dXdv12=0;
346     dXdv22=0;
347 gross 776 for (s=0;s<numShape; s++) {
348 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
349     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
350     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
351     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
352     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
353     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
354     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
355     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
356     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
357     dXdv02+=X0_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
358     dXdv12+=X1_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
359     dXdv22+=X2_loc*DSDv[INDEX3(s,2,q,numShape,LOCDIM)];
360 gross 776 }
361 gross 777 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
362 gross 776 if (D==0.) {
363 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
364 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
365 gross 776 } else {
366     invD=1./D;
367 gross 777 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
368     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
369     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
370     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
371     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
372     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
373     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
374     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
375     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
376    
377 gross 776 for (s=0;s<numTest; s++) {
378 gross 781 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=
379     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX20;
380     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=
381     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX21;
382     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=
383     DTDv[INDEX3(s,0,q,numTest,LOCDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numTest,LOCDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numTest,LOCDIM)]*dvdX22;
384 gross 776 }
385     }
386 gross 777 m0=dXdv10*dXdv21-dXdv20*dXdv11;
387     m1=dXdv20*dXdv01-dXdv00*dXdv21;
388     m2=dXdv00*dXdv11-dXdv10*dXdv01;
389     volume[INDEX2(q,e,numQuad)]=sqrt(m0*m0+m1*m1+m2*m2)*QuadWeights[q];
390 gross 776 }
391 jgs 82 }
392    
393 gross 776 }
394     #undef DIM
395 gross 781 #undef LOCDIM
396 gross 776 }
397     /**************************************************************/
398     /* */
399 gross 777 /* Jacobean 2D manifold in 3D with 2D elements */
400 gross 776 /* */
401 gross 777 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad,double* QuadWeights,
402     dim_t numShape, dim_t numElements, dim_t numNodes, index_t* nodes,
403     double* DSDv, dim_t numTest,double* DTDv,
404     double* dTdX, double* volume, index_t* element_id) {
405     #define DIM 3
406 gross 781 #define LOCDIM 2
407 gross 776 register int e,q,s;
408     char error_msg[LenErrorMsg_MAX];
409     #pragma omp parallel
410 jgs 82 {
411 gross 777 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
412 gross 853 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
413     X0_loc, X1_loc, X2_loc;
414 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)
415 gross 776 for(e=0;e<numElements;e++){
416     for (q=0;q<numQuad;q++) {
417     dXdv00=0;
418     dXdv10=0;
419 gross 777 dXdv20=0;
420     dXdv01=0;
421     dXdv11=0;
422     dXdv21=0;
423 gross 776 for (s=0;s<numShape; s++) {
424 gross 853 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
425     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
426     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
427     dXdv00+=X0_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
428     dXdv10+=X1_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
429     dXdv20+=X2_loc*DSDv[INDEX3(s,0,q,numShape,LOCDIM)];
430     dXdv01+=X0_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
431     dXdv11+=X1_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
432     dXdv21+=X2_loc*DSDv[INDEX3(s,1,q,numShape,LOCDIM)];
433 gross 776 }
434 gross 777 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
435     m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
436     m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
437     D=m00*m11-m01*m01;
438 gross 776 if (D==0.) {
439 gross 777 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
440 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
441 gross 776 } else {
442     invD=1./D;
443 gross 777 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
444     dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
445     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
446     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
447     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
448     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
449 gross 776 for (s=0;s<numTest; s++) {
450 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;
451     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;
452     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;
453 gross 776 }
454     volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
455     }
456     }
457     }
458 jgs 82
459     }
460 gross 776 #undef DIM
461 gross 781 #undef LOCDIM
462 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26