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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1628 - (hide annotations)
Fri Jul 11 13:12:46 2008 UTC (10 years, 11 months ago) by phornby
Original Path: trunk/finley/src/Assemble_jacobeans.c
File MIME type: text/plain
File size: 44283 byte(s)

Merge in /branches/windows_from_1456_trunk_1620_merged_in branch.

You will find a preserved pre-merge trunk in tags under tags/trunk_at_1625.
That will be useful for diffing & checking on my stupidity.

Here is a list of the conflicts and their resolution at this
point in time.


=================================================================================
(LLWS == looks like white space).

finley/src/Assemble_addToSystemMatrix.c - resolve to branch - unused var. may be wrong.....
finley/src/CPPAdapter/SystemMatrixAdapter.cpp - resolve to branch - LLWS
finley/src/CPPAdapter/MeshAdapter.cpp - resolve to branch - LLWS
paso/src/PCG.c - resolve to branch - unused var fixes.
paso/src/SolverFCT.c - resolve to branch - LLWS
paso/src/FGMRES.c - resolve to branch - LLWS
paso/src/Common.h - resolve to trunk version. It's omp.h's include... not sure it's needed,
but for the sake of saftey.....
paso/src/Functions.c - resolve to branch version, indentation/tab removal and return error
on bad unimplemented Paso_FunctionCall.
paso/src/SolverFCT_solve.c - resolve to branch version, unused vars
paso/src/SparseMatrix_MatrixVector.c - resolve to branch version, unused vars.
escript/src/Utils.cpp - resloved to branch, needs WinSock2.h
escript/src/DataExpanded.cpp - resolved to branch version - LLWS
escript/src/DataFactory.cpp - resolve to branch version
=================================================================================

This currently passes tests on linux (debian), but is not checked on windows or Altix yet.

This checkin is to make a trunk I can check out for windows to do tests on it.

Known outstanding problem is in the operator=() method of exceptions
causing warning messages on the intel compilers.

May the God of doughnuts have mercy on my soul.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26