/[escript]/branches/domexper/dudley/src/Assemble_jacobeans.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3184 - (show annotations)
Wed Sep 15 00:23:42 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 11793 byte(s)
Removed volume from ElementFile_Jacobeans

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * 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
14
15 #include "Assemble.h"
16 #include "Util.h"
17 #ifdef _OPENMP
18 #include <omp.h>
19 #endif
20
21
22 /* input:
23
24 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 #include "ShapeTable.h"
44
45 #define SCALING(_nsub_,_dim_) pow(1./(double)(_nsub_),1./(double)(_dim_))
46
47 /**************************************************************/
48 /* */
49 /* Jacobean 2D with area element */
50 /* */
51 void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
52 double* dTdX, double* absD, double* quadweight, index_t* element_id)
53 {
54 #define DIM 2
55 #define LOCDIM 2
56 register int e,q,s;
57 char error_msg[LenErrorMsg_MAX];
58 *quadweight=(numQuad==1)?1./2:1./6; /* numQuad is 1 or 3 */
59 const dim_t numTest=3; // hoping this is used in constant folding
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 #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 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 absD[e]=ABS(D);
85 if (D==0.)
86 {
87 sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
88 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
89 }
90 else
91 {
92 invD=1./D;
93 dvdX00= dXdv11*invD;
94 dvdX10=-dXdv10*invD;
95 dvdX01=-dXdv01*invD;
96 dvdX11= dXdv00*invD;
97 if (numQuad==1)
98 {
99 for (s=0;s<3; s++)
100 {
101 #define DTDXSET(P,Q) dTdX[INDEX4(s,P,Q,e,numTest,DIM,numQuad)] = DTDV_2D[s][0]*dvdX0##P+DTDV_2D[s][1]*dvdX1##P
102
103 DTDXSET(0,0);
104 DTDXSET(1,0);
105 }
106 // volume[INDEX2(0,e,1)]=ABS(D)*(*quadweight);
107 }
108 else // numQuad==3
109 {
110 for (q=0;q<numTest;++q) // relying on unroll loops to optimise this
111 {
112 for (s=0;s<3; s++)
113 {
114 DTDXSET(0,q);
115 DTDXSET(1,q);
116 }
117 // volume[INDEX2(q,e,3)]=ABS(D)*(*quadweight);
118 }
119 }
120 }
121 }
122 } // end parallel
123 #undef DIM
124 #undef LOCDIM
125 #undef DTDXSET
126 #undef COMPDXDV0
127 #undef COMPDXDV1
128 }
129 /**************************************************************/
130 /* */
131 /* Jacobean 1D manifold in 2D and 1D elements */
132 /* */
133 void Assemble_jacobeans_2D_M1D_E1D(double* coordinates, dim_t numQuad,
134 dim_t numElements, dim_t numNodes, index_t* nodes,
135 double* dTdX, double* absD, double* quadweight, index_t* element_id)
136 {
137 #define DIM 2
138 #define LOCDIM 1
139 register int e;
140 char error_msg[LenErrorMsg_MAX];
141 const dim_t numTest=2;
142 *quadweight=(numQuad==1)?1.0:0.5;
143 // numQuad is 1 or 2
144 #pragma omp parallel
145 {
146 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
147 #pragma omp for private(e,q,s,dXdv00,dXdv10,dvdX00,dvdX01,D,invD,X0_loc, X1_loc) schedule(static)
148 for(e=0;e<numElements;e++)
149 {
150 dXdv00=0;
151 dXdv10=0;
152 dXdv00+=coordinates[INDEX2(0,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(0,nodes[INDEX2(1,e,numNodes)],DIM)];
153 dXdv00+=coordinates[INDEX2(1,nodes[INDEX2(0,e,numNodes)],DIM)]*(-1.)+coordinates[INDEX2(1,nodes[INDEX2(1,e,numNodes)],DIM)];
154 D=dXdv00*dXdv00+dXdv10*dXdv10;
155 if (D==0.)
156 {
157 sprintf(error_msg,"Assemble_jacobeans_2D_M1D_E1D: element %d (id %d) has length zero.",e,element_id[e]);
158 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
159 } else {
160 invD=1./D;
161 dvdX00=dXdv00*invD;
162 dvdX01=dXdv10*invD;
163 // The number of quad points is 1 or 2
164 dTdX[INDEX4(0,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
165 dTdX[INDEX4(0,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
166 dTdX[INDEX4(1,0,0,e,numTest,DIM,numQuad)]=-1*dvdX00;
167 dTdX[INDEX4(1,1,0,e,numTest,DIM,numQuad)]=-1*dvdX01;
168 // volume[INDEX2(0,e,numQuad)]=sqrt(D)*(*quadweight);
169 absD[e]=sqrt(D);
170 if (numQuad==2)
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 }
180 } // end parallel
181 #undef DIM
182 #undef LOCDIM
183 }
184
185 /**************************************************************/
186 /* */
187 /* Jacobean 3D */
188 /* */
189 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
190 double* dTdX, double* absD, double* quadweight, index_t* element_id)
191 {
192 #define DIM 3
193 #define LOCDIM 3
194 int e,q,s;
195 char error_msg[LenErrorMsg_MAX];
196 // numQuad is 1 or 4
197 const dim_t numShape=4, numTest=4;
198 *quadweight=(numQuad==1)?1./6:1./24;
199
200 #pragma omp parallel
201 {
202 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
203 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD,
204 X0_loc,X1_loc,X2_loc;
205 #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)
206 for(e=0;e<numElements;e++)
207 {
208 dXdv00=0;
209 dXdv10=0;
210 dXdv20=0;
211 dXdv01=0;
212 dXdv11=0;
213 dXdv21=0;
214 dXdv02=0;
215 dXdv12=0;
216 dXdv22=0;
217 for (s=0;s<numShape; s++)
218 {
219 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
220 X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
221 X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
222 dXdv00+=X0_loc*DTDV_3D[s][0];
223 dXdv10+=X1_loc*DTDV_3D[s][0];
224 dXdv20+=X2_loc*DTDV_3D[s][0];
225 dXdv01+=X0_loc*DTDV_3D[s][1];
226 dXdv11+=X1_loc*DTDV_3D[s][1];
227 dXdv21+=X2_loc*DTDV_3D[s][1];
228 dXdv02+=X0_loc*DTDV_3D[s][2];
229 dXdv12+=X1_loc*DTDV_3D[s][2];
230 dXdv22+=X2_loc*DTDV_3D[s][2];
231 }
232 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
233 absD[e]=ABS(D);
234 if (D==0.)
235 {
236 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
237 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
238 } else {
239 invD=1./D;
240 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
241 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
242 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
243 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
244 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
245 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
246 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
247 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
248 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
249 for (q=0;q<numQuad;q++)
250 {
251 for (s=0;s<numTest; s++)
252 {
253 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV_3D[s][0]*dvdX00+DTDV_3D[s][1]*dvdX10+DTDV_3D[s][2]*dvdX20;
254 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV_3D[s][0]*dvdX01+DTDV_3D[s][1]*dvdX11+DTDV_3D[s][2]*dvdX21;
255 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV_3D[s][0]*dvdX02+DTDV_3D[s][1]*dvdX12+DTDV_3D[s][2]*dvdX22;
256 }
257 // volume[INDEX2(q,e,numQuad)]=ABS(D)*(*quadweight);
258 }
259 }
260 }
261 } // end parallel
262 #undef DIM
263 #undef LOCDIM
264 }
265
266 /**************************************************************/
267 /* */
268 /* Jacobean 2D manifold in 3D with 2D elements */
269 /* */
270 void Assemble_jacobeans_3D_M2D_E2D(double* coordinates, dim_t numQuad, dim_t numElements, dim_t numNodes, index_t* nodes,
271 double* dTdX, double* absD, double* quadweight, index_t* element_id)
272 {
273 #define DIM 3
274 #define LOCDIM 2
275 register int e,q,s;
276 char error_msg[LenErrorMsg_MAX];
277 // numQuad is 1 or 3
278 *quadweight=(numQuad==1)?1./2:1./6;
279 const double DTDV[3][2]={{-1.,-1.},{1.,0.},{0.,1.}};
280 const dim_t numShape=3, numTest=3;
281 #pragma omp parallel
282 {
283 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
284 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD,
285 X0_loc, X1_loc, X2_loc;
286 #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)
287 for(e=0;e<numElements;e++)
288 {
289 dXdv00=0;
290 dXdv10=0;
291 dXdv20=0;
292 dXdv01=0;
293 dXdv11=0;
294 dXdv21=0;
295 for (s=0;s<numShape; s++)
296 {
297 X0_loc=coordinates[INDEX2(0,nodes[INDEX2(s,e,numNodes)],DIM)];
298 X1_loc=coordinates[INDEX2(1,nodes[INDEX2(s,e,numNodes)],DIM)];
299 X2_loc=coordinates[INDEX2(2,nodes[INDEX2(s,e,numNodes)],DIM)];
300 dXdv00+=X0_loc*DTDV[s][0];
301 dXdv10+=X1_loc*DTDV[s][0];
302 dXdv20+=X2_loc*DTDV[s][0];
303 dXdv01+=X0_loc*DTDV[s][1];
304 dXdv11+=X1_loc*DTDV[s][1];
305 dXdv21+=X2_loc*DTDV[s][1];
306 }
307 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
308 m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
309 m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
310 D=m00*m11-m01*m01;
311 absD[e]=sqrt(D);
312 if (D==0.)
313 {
314 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
315 Dudley_setError(ZERO_DIVISION_ERROR,error_msg);
316 } else {
317 invD=1./D;
318 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
319 dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
320 dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
321 dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
322 dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
323 dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
324 for (q=0;q<numQuad;q++)
325 {
326 for (s=0;s<numTest; s++)
327 {
328 dTdX[INDEX4(s,0,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX00+DTDV[s][1]*dvdX10;
329 dTdX[INDEX4(s,1,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX01+DTDV[s][1]*dvdX11;
330 dTdX[INDEX4(s,2,q,e,numTest,DIM,numQuad)]=DTDV[s][0]*dvdX02+DTDV[s][1]*dvdX12;
331 }
332 // volume[INDEX2(q,e,numQuad)]=sqrt(D)*(*quadweight);
333 }
334 }
335 }
336 } // end parallel section
337 #undef DIM
338 #undef LOCDIM
339 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26