/[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 3173 - (hide annotations)
Fri Sep 10 03:24:38 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 11782 byte(s)
Indented Assemble_jacobeans consistantly.
Completed the first pass on the remainder of the functions in Assemble_jacobeans.

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 jfenwick 3173 }
120 jfenwick 3169 }
121 jfenwick 3173 } // end parallel
122 gross 776 #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 3172 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,
133     dim_t numElements, dim_t numNodes, index_t* nodes,
134 jfenwick 3173 double* dTdX, double* volume, index_t* element_id)
135     {
136     #define DIM 2
137     #define LOCDIM 1
138     register int e;
139     char error_msg[LenErrorMsg_MAX];
140     const dim_t numTest=2;
141     double quadweight=(numQuad==1)?1.0:0.5;
142     // numQuad is 1 or 2
143     #pragma omp parallel
144     {
145     register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
146     #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
147     for(e=0;e<numElements;e++)
148     {
149     dXdv00=0;
150     dXdv10=0;
151     dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
152     dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
153     D=dXdv00*dXdv00+dXdv10*dXdv10;
154     if (D==0.)
155     {
156     sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
157     Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
158     } else {
159     invD=1./D;
160     dvdX00=dXdv00*invD;
161     dvdX01=dXdv10*invD;
162     // The number of quad points is 1 or 2
163     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 jfenwick 3171 if (numQuad==2)
169     {
170 jfenwick 3173 dTdX[INDEX4(0,0,1,e,numTest,DIM,numQuad)]=dvdX00;
171     dTdX[INDEX4(0,1,1,e,numTest,DIM,numQuad)]=dvdX01;
172     dTdX[INDEX4(1,0,1,e,numTest,DIM,numQuad)]=dvdX00;
173     dTdX[INDEX4(1,1,1,e,numTest,DIM,numQuad)]=dvdX01;
174     volume[INDEX2(1,e,numQuad)]=sqrt(D)*quadweight;
175 jfenwick 3171 }
176 jfenwick 3173 }
177     }
178     } // end parallel
179     #undef DIM
180     #undef LOCDIM
181 gross 781 }
182 jfenwick 3171
183 gross 781 /**************************************************************/
184     /* */
185 gross 777 /* Jacobean 3D */
186 gross 776 /* */
187 jfenwick 3172 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
188 jfenwick 3173 double* dTdX, double* volume, index_t* element_id)
189     {
190     #define DIM 3
191     #define LOCDIM 3
192     int e,q,s;
193     char error_msg[LenErrorMsg_MAX];
194     // numQuad is 1 or 4
195     const dim_t numShape=4, numTest=4;
196     const double quadweight=(numQuad==1)?1./6:1./24;
197     const double DTDV[4][3]={{-1, -1, -1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
198     #pragma omp parallel
199     {
200     register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
201 gross 853 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
202     X0_loc,X1_loc,X2_loc;
203 jfenwick 3173 #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)
204     for(e=0;e<numElements;e++)
205     {
206     dXdv00=0;
207     dXdv10=0;
208     dXdv20=0;
209     dXdv01=0;
210     dXdv11=0;
211     dXdv21=0;
212     dXdv02=0;
213     dXdv12=0;
214     dXdv22=0;
215     for (s=0;s<numShape; s++)
216     {
217     X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
218     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
219     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
220     dXdv00+=X0_loc*DTDV[s][0];
221     dXdv10+=X1_loc*DTDV[s][0];
222     dXdv20+=X2_loc*DTDV[s][0];
223     dXdv01+=X0_loc*DTDV[s][1];
224     dXdv11+=X1_loc*DTDV[s][1];
225     dXdv21+=X2_loc*DTDV[s][1];
226     dXdv02+=X0_loc*DTDV[s][2];
227     dXdv12+=X1_loc*DTDV[s][2];
228     dXdv22+=X2_loc*DTDV[s][2];
229     }
230     D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
231     if (D==0.)
232     {
233     sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
234     Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
235     } else {
236     invD=1./D;
237     dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
238     dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
239     dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
240     dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
241     dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
242     dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
243     dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
244     dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
245     dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
246     for (q=0;q<numQuad;q++)
247     {
248     for (s=0;s<numTest; s++)
249     {
250     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX00+DTDV[s][1]*dvdX10+DTDV[s][2]*dvdX20;
251     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX01+DTDV[s][1]*dvdX11+DTDV[s][2]*dvdX21;
252     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX02+DTDV[s][1]*dvdX12+DTDV[s][2]*dvdX22;
253     }
254     volume[INDEX2(q,e,numQuad)]=ABS(D)*quadweight;
255     }
256     }
257     }
258     } // end parallel
259 gross 776 #undef DIM
260 gross 781 #undef LOCDIM
261 gross 776 }
262 gross 777
263 gross 776 /**************************************************************/
264     /* */
265 gross 777 /* Jacobean 2D manifold in 3D with 2D elements */
266 gross 776 /* */
267 jfenwick 3173 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
268     double* dTdX, double* volume, index_t* element_id)
269     {
270     #define DIM 3
271     #define LOCDIM 2
272     register int e,q,s;
273     char error_msg[LenErrorMsg_MAX];
274     // numQuad is 1 or 3
275     const double quadweight=(numQuad==1)?1./2:1./6;
276     const double DTDV[3][2]={{-1.,-1.},{1.,0.},{0.,1.}};
277     const dim_t numShape=3, numTest=3;
278     #pragma omp parallel
279     {
280     register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
281 gross 853 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
282     X0_loc, X1_loc, X2_loc;
283 jfenwick 3173 #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)
284     for(e=0;e<numElements;e++)
285     {
286     dXdv00=0;
287     dXdv10=0;
288     dXdv20=0;
289     dXdv01=0;
290     dXdv11=0;
291     dXdv21=0;
292     for (s=0;s<numShape; s++)
293     {
294     X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
295     X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
296     X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
297     dXdv00+=X0_loc*DTDV[s][0];
298     dXdv10+=X1_loc*DTDV[s][0];
299     dXdv20+=X2_loc*DTDV[s][0];
300     dXdv01+=X0_loc*DTDV[s][1];
301     dXdv11+=X1_loc*DTDV[s][1];
302     dXdv21+=X2_loc*DTDV[s][1];
303     }
304     m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
305     m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
306     m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
307     D=m00*m11-m01*m01;
308     if (D==0.)
309     {
310     sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
311     Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
312     } else {
313     invD=1./D;
314     dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
315     dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
316     dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
317     dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
318     dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
319     dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
320     for (q=0;q<numQuad;q++)
321     {
322     for (s=0;s<numTest; s++)
323     {
324     dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX00+DTDV[s][1]*dvdX10;
325     dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX01+DTDV[s][1]*dvdX11;
326     dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX02+DTDV[s][1]*dvdX12;
327     }
328     volume[INDEX2(q,e,numQuad)]=sqrt(D)*quadweight;
329     }
330     }
331     }
332     } // end parallel section
333     #undef DIM
334     #undef LOCDIM
335 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26