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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26