/[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 3170 - (show annotations)
Thu Sep 9 05:40:25 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 44821 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26