/[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 3169 - (show annotations)
Thu Sep 9 03:21:35 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 45750 byte(s)
Some moving

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26