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

Contents of /temp/finley/src/Assemble_jacobeans.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 776 - (show annotations)
Wed Jul 12 00:07:31 2006 UTC (13 years, 3 months ago) by gross
Original Path: trunk/finley/src/Assemble_jacobeans.c
File MIME type: text/plain
File size: 17351 byte(s)
basic code for persistence of jacobeans added. routines are not called yet
1 /*
2 ************************************************************
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 */
12
13
14 /**************************************************************/
15
16 /* Author: gross@access.edu.au */
17 /* Version: $Id$ */
18
19 /**************************************************************/
20
21 #include "Assemble.h"
22 #include "Util.h"
23 #ifdef _OPENMP
24 #include <omp.h>
25 #endif
26
27 /**************************************************************/
28 /* */
29 /* Jacobean 1D */
30 /* */
31 void Assemble_jacobeans_1D(double* coordinates, dim_t numQuad,double* QuadWeights,
32 dim_t numShape, dim_t numElements,index_t* nodes,
33 double* DSDv, dim_t numTest,double* DTDv,
34 double* dTdX, double* volume, index_t* element_id) {
35 #define DIM 1
36 #define LOCALDIM 1
37 register int e,q,s;
38 char error_msg[LenErrorMsg_MAX];
39 #pragma omp parallel
40 {
41 register double D,invD;
42 double X0[numShape];
43 #pragma omp for private(e,q,s) schedule(static)
44 for(e=0;e<numElements;e++){
45 for (s=0;s<numShape; s++) X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
46 for (q=0;q<numQuad;q++) {
47 D=0;
48 for (s=0;s<numShape; s++) D+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
49 if (D==0.) {
50 sprintf(error_msg,"Assemble_jacobeans_1D: element %d (id %d) has length zero.",e,element_id[e]);
51 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
52 } else {
53 invD=1./D;
54 for (s=0;s<numTest; s++) dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*invD;
55 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
56 }
57 }
58 }
59
60 }
61 #undef DIM
62 #undef LOCALDIM
63
64 }
65 /**************************************************************/
66 /* */
67 /* Jacobean 2D */
68 /* */
69 void Assemble_jacobeans_2D(double* coordinates, dim_t numQuad,double* QuadWeights,
70 dim_t numShape, dim_t numElements,index_t* nodes,
71 double* DSDv, dim_t numTest,double* DTDv,
72 double* dTdX, double* volume, index_t* element_id) {
73 #define DIM 2
74 #define LOCALDIM 2
75 register int e,q,s;
76 char error_msg[LenErrorMsg_MAX];
77 #pragma omp parallel
78 {
79 register double dXdv00,dXdv10,dXdv01,dXdv11,
80 dvdX00,dvdX10,dvdX01,dvdX11, D,invD;
81 double X0[numShape], X1[numShape];
82 #pragma omp for private(e,q,s) schedule(static)
83 for(e=0;e<numElements;e++){
84 for (s=0;s<numShape; s++) {
85 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
86 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
87 }
88 for (q=0;q<numQuad;q++) {
89 dXdv00=0;
90 dXdv10=0;
91 dXdv01=0;
92 dXdv11=0;
93 for (s=0;s<numShape; s++) {
94 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
95 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
96 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
97 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
98 }
99 D = dXdv00*dXdv11 - dXdv01*dXdv10;
100 if (D==0.) {
101 sprintf(error_msg,"Assemble_jacobeans_2D: element %d (id %d) has area zero.",e,element_id[e]);
102 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
103 } else {
104 invD=1./D;
105 dvdX00= dXdv11*invD;
106 dvdX10=-dXdv10*invD;
107 dvdX01=-dXdv01*invD;
108 dvdX11= dXdv00*invD;
109
110 for (s=0;s<numTest; s++) {
111 dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
112 dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
113 }
114 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
115 }
116 }
117 }
118
119 }
120 #undef DIM
121 #undef LOCALDIM
122 }
123 /**************************************************************/
124 /* */
125 /* Jacobean 3D */
126 /* */
127 void Assemble_jacobeans_3D(double* coordinates, dim_t numQuad,double* QuadWeights,
128 dim_t numShape, dim_t numElements,index_t* nodes,
129 double* DSDv, dim_t numTest,double* DTDv,
130 double* dTdX, double* volume, index_t* element_id) {
131 #define DIM 3
132 #define LOCALDIM 3
133 register int e,q,s;
134 char error_msg[LenErrorMsg_MAX];
135 #pragma omp parallel
136 {
137 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,dXdv02,dXdv12,dXdv22,
138 dvdX00,dvdX10,dvdX20,dvdX01,dvdX11,dvdX21,dvdX02,dvdX12,dvdX22, D,invD;
139 double X0[numShape], X1[numShape], X2[numShape];
140 #pragma omp for private(e,q,s) schedule(static)
141 for(e=0;e<numElements;e++){
142 for (s=0;s<numShape; s++) {
143 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
144 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
145 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
146 }
147 for (q=0;q<numQuad;q++) {
148 dXdv00=0;
149 dXdv10=0;
150 dXdv20=0;
151 dXdv01=0;
152 dXdv11=0;
153 dXdv21=0;
154 dXdv02=0;
155 dXdv12=0;
156 dXdv22=0;
157 for (s=0;s<numShape; s++) {
158 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
159 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
160 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
161 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
162 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
163 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
164 dXdv02+=X0[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
165 dXdv12+=X1[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
166 dXdv22+=X2[s]*DSDv[INDEX3(s,2,q,numShape,LOCALDIM)];
167 }
168 D = dXdv00*(dXdv11*dXdv22-dXdv12*dXdv21)+ dXdv01*(dXdv20*dXdv12-dXdv10*dXdv22)+dXdv02*(dXdv10*dXdv21-dXdv20*dXdv11);
169 if (D==0.) {
170 sprintf(error_msg,"Assemble_jacobeans_3D: element %d (id %d) has volume zero.",e,element_id[e]);
171 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
172 } else {
173 invD=1./D;
174 dvdX00=(dXdv11*dXdv22-dXdv12*dXdv21)*invD;
175 dvdX10=(dXdv20*dXdv12-dXdv10*dXdv22)*invD;
176 dvdX20=(dXdv10*dXdv21-dXdv20*dXdv11)*invD;
177 dvdX01=(dXdv02*dXdv21-dXdv01*dXdv22)*invD;
178 dvdX11=(dXdv00*dXdv22-dXdv20*dXdv02)*invD;
179 dvdX21=(dXdv01*dXdv20-dXdv00*dXdv21)*invD;
180 dvdX02=(dXdv01*dXdv12-dXdv02*dXdv11)*invD;
181 dvdX12=(dXdv02*dXdv10-dXdv00*dXdv12)*invD;
182 dvdX22=(dXdv00*dXdv11-dXdv01*dXdv10)*invD;
183
184 for (s=0;s<numTest; s++) {
185 dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX20;
186 dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX21;
187 dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12+DTDv[INDEX3(s,2,q,numShape,LOCALDIM)]*dvdX22;
188 }
189 volume[INDEX2(q,e,numQuad)]=ABS(D)*QuadWeights[q];
190 }
191 }
192 }
193
194 }
195 #undef DIM
196 #undef LOCALDIM
197 }
198 /**************************************************************/
199 /* */
200 /* Jacobean 2D manifold in 3D */
201 /* */
202 void Assemble_jacobeans_3D_M2D(double* coordinates, dim_t numQuad,double* QuadWeights,
203 dim_t numShape, dim_t numElements,index_t* nodes,
204 double* DSDv, dim_t numTest,double* DTDv,
205 double* dTdX, double* volume, index_t* element_id) {
206 #define DIM 3
207 #define LOCALDIM 2
208 register int e,q,s;
209 char error_msg[LenErrorMsg_MAX];
210 #pragma omp parallel
211 {
212 register double dXdv00,dXdv10,dXdv20,dXdv01,dXdv11,dXdv21,m00,m01,m11,
213 dvdX00,dvdX01,dvdX02,dvdX10,dvdX11,dvdX12,D,invD;
214 double X0[numShape], X1[numShape], X2[numShape];
215 #pragma omp for private(e,q,s) schedule(static)
216 for(e=0;e<numElements;e++){
217 for (s=0;s<numShape; s++) {
218 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
219 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
220 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
221 }
222 for (q=0;q<numQuad;q++) {
223 dXdv00=0;
224 dXdv10=0;
225 dXdv20=0;
226 dXdv01=0;
227 dXdv11=0;
228 dXdv21=0;
229 for (s=0;s<numShape; s++) {
230 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
231 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
232 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
233 dXdv01+=X0[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
234 dXdv11+=X1[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
235 dXdv21+=X2[s]*DSDv[INDEX3(s,1,q,numShape,LOCALDIM)];
236 }
237 m00=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
238 m01=dXdv00*dXdv01+dXdv10*dXdv11+dXdv20*dXdv21;
239 m11=dXdv01*dXdv01+dXdv11*dXdv11+dXdv21*dXdv21;
240 D=m00*m11-m01*m01;
241 if (D==0.) {
242 sprintf(error_msg,"Assemble_jacobeans_3D_M2D: element %d (id %d) has area zero.",e,element_id[e]);
243 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
244 } else {
245 invD=1./D;
246 dvdX00=( m00*dXdv00-m01*dXdv01)*invD;
247 dvdX01=( m00*dXdv10-m01*dXdv11)*invD;
248 dvdX02=( m00*dXdv20-m01*dXdv21)*invD;
249 dvdX10=(-m01*dXdv00+m11*dXdv01)*invD;
250 dvdX11=(-m01*dXdv10+m11*dXdv11)*invD;
251 dvdX12=(-m01*dXdv20+m11*dXdv21)*invD;
252 for (s=0;s<numTest; s++) {
253 dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX10;
254 dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX11;
255 dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02+DTDv[INDEX3(s,1,q,numShape,LOCALDIM)]*dvdX12;
256 }
257 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
258 }
259 }
260 }
261
262 }
263 #undef DIM
264 #undef LOCALDIM
265 }
266 /**************************************************************/
267 /* */
268 /* Jacobean 1D manifold in 3D */
269 /* */
270 void Assemble_jacobeans_3D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,
271 dim_t numShape, dim_t numElements,index_t* nodes,
272 double* DSDv, dim_t numTest,double* DTDv,
273 double* dTdX, double* volume, index_t* element_id) {
274 #define DIM 3
275 #define LOCALDIM 1
276 register int e,q,s;
277 char error_msg[LenErrorMsg_MAX];
278 #pragma omp parallel
279 {
280 register double dXdv00,dXdv10,dXdv20,dvdX00,dvdX01,dvdX02,D,invD;
281 double X0[numShape], X1[numShape], X2[numShape];
282 #pragma omp for private(e,q,s) schedule(static)
283 for(e=0;e<numElements;e++){
284 for (s=0;s<numShape; s++) {
285 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
286 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
287 X2[s]=coordinates[INDEX2(2,nodes[INDEX2(s,e,numShape)],DIM)];
288 }
289 for (q=0;q<numQuad;q++) {
290 dXdv00=0;
291 dXdv10=0;
292 dXdv20=0;
293 for (s=0;s<numShape; s++) {
294 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
295 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
296 dXdv20+=X2[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
297 }
298 D=dXdv00*dXdv00+dXdv10*dXdv10+dXdv20*dXdv20;
299 if (D==0.) {
300 sprintf(error_msg,"Assemble_jacobeans_3D_M1D: element %d (id %d) has length zero.",e,element_id[e]);
301 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
302 } else {
303 invD=1./D;
304 dvdX00=dXdv00*invD;
305 dvdX01=dXdv10*invD;
306 dvdX02=dXdv20*invD;
307 for (s=0;s<numTest; s++) {
308 dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;
309 dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;
310 dTdX[INDEX4(s,2,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX02;
311 }
312 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
313 }
314 }
315 }
316
317 }
318 #undef DIM
319 #undef LOCALDIM
320 }
321 /**************************************************************/
322 /* */
323 /* Jacobean 1D manifold in 2D */
324 /* */
325 void Assemble_jacobeans_2D_M1D(double* coordinates, dim_t numQuad,double* QuadWeights,
326 dim_t numShape, dim_t numElements,index_t* nodes,
327 double* DSDv, dim_t numTest,double* DTDv,
328 double* dTdX, double* volume, index_t* element_id) {
329 #define DIM 2
330 #define LOCALDIM 1
331 register int e,q,s;
332 char error_msg[LenErrorMsg_MAX];
333 #pragma omp parallel
334 {
335 register double dXdv00,dXdv10,dvdX00,dvdX01,D,invD;
336 double X0[numShape], X1[numShape];
337 #pragma omp for private(e,q,s) schedule(static)
338 for(e=0;e<numElements;e++){
339 for (s=0;s<numShape; s++) {
340 X0[s]=coordinates[INDEX2(0,nodes[INDEX2(s,e,numShape)],DIM)];
341 X1[s]=coordinates[INDEX2(1,nodes[INDEX2(s,e,numShape)],DIM)];
342 }
343 for (q=0;q<numQuad;q++) {
344 dXdv00=0;
345 dXdv10=0;
346 for (s=0;s<numShape; s++) {
347 dXdv00+=X0[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
348 dXdv10+=X1[s]*DSDv[INDEX3(s,0,q,numShape,LOCALDIM)];
349 }
350 D=dXdv00*dXdv00+dXdv10*dXdv10;
351 if (D==0.) {
352 sprintf(error_msg,"Assemble_jacobeans_2D_M1D: element %d (id %d) has length zero.",e,element_id[e]);
353 Finley_setError(ZERO_DIVISION_ERROR,error_msg);
354 } else {
355 invD=1./D;
356 dvdX00=dXdv00*invD;
357 dvdX01=dXdv10*invD;
358 for (s=0;s<numTest; s++) {
359 dTdX[INDEX4(s,0,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX00;
360 dTdX[INDEX4(s,1,q,e,numTest,LOCALDIM,numQuad)]=DTDv[INDEX3(s,0,q,numShape,LOCALDIM)]*dvdX01;
361 }
362 volume[INDEX2(q,e,numQuad)]=sqrt(D)*QuadWeights[q];
363 }
364 }
365 }
366
367 }
368 #undef DIM
369 #undef LOCALDIM
370 }
371 /*
372 * $Log$
373 *
374 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26