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

Annotation of /trunk/finley/src/Assemble_gradient.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1028 - (hide annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 10 months ago) by gross
File MIME type: text/plain
File size: 18121 byte(s)
modifications to be compliant with _WIN32. The substitutes for asinh, acosh, atanh are still missing (erf will through an exception)
1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 82
13     /**************************************************************/
14    
15 gross 781 /* assemblage rjacines: calculate the gradient of nodal data at quadrature points */
16 jgs 82
17     /**************************************************************/
18    
19 jgs 150 /* Copyrights by ACcESS Australia, 2003,2004,2005 */
20 jgs 82 /* author: gross@access.edu.au */
21 jgs 150 /* version: $Id$ */
22 jgs 82
23     /**************************************************************/
24    
25     #include "Assemble.h"
26     #include "Util.h"
27     #ifdef _OPENMP
28     #include <omp.h>
29     #endif
30     /*****************************************************************/
31    
32    
33     void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements,
34 gross 781 escriptDataC* grad_data,escriptDataC* data) {
35 jgs 82
36 gross 1028 dim_t numNodes, numShapes, numLocalNodes, numComps, NN;
37 jgs 123 type_t data_type=getFunctionSpaceType(data);
38 gross 781 type_t grad_data_type=getFunctionSpaceType(grad_data);
39     bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
40 gross 1028 index_t dof_offset, s_offset;
41     Finley_ElementFile_Jacobeans* jac=NULL;
42     Finley_resetError();
43     if (nodes==NULL || elements==NULL) return;
44     numComps=getDataPointSize(data);
45     NN=elements->ReferenceElement->Type->numNodes;
46 jgs 82
47     if (data_type==FINLEY_NODES) {
48 gross 781 reducedShapefunction=FALSE;
49 jgs 82 numNodes=nodes->numNodes;
50 gross 781 numShapes=elements->ReferenceElement->Type->numShapes;
51     numLocalNodes=elements->ReferenceElement->Type->numNodes;
52 bcumming 751 }
53 gross 781 /* lock these two options jac for the MPI version */
54     #ifndef PASO_MPI
55     else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
56     reducedShapefunction=FALSE;
57 jgs 82 numNodes=nodes->numDegreesOfFreedom;
58 gross 781 numShapes=elements->ReferenceElement->Type->numShapes;
59     numLocalNodes=elements->ReferenceElement->Type->numNodes;
60 jgs 82 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
61 gross 781 reducedShapefunction=TRUE;
62 jgs 82 numNodes=nodes->reducedNumDegreesOfFreedom;
63 gross 781 numShapes=elements->LinearReferenceElement->Type->numShapes;
64     numLocalNodes=elements->LinearReferenceElement->Type->numNodes;
65     }
66     #endif
67     else {
68     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
69 jgs 82 }
70 gross 781 if (grad_data_type==FINLEY_ELEMENTS) {
71     reducedIntegrationOrder=FALSE;
72     dof_offset=0;
73     } else if (grad_data_type==FINLEY_FACE_ELEMENTS) {
74     reducedIntegrationOrder=FALSE;
75     dof_offset=0;
76     } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_1) {
77     reducedIntegrationOrder=FALSE;
78     dof_offset=0;
79     } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2) {
80     reducedIntegrationOrder=FALSE;
81     /* don't reset offset */
82 jgs 82 } else {
83 gross 781 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculated for requested location.");
84 jgs 82 }
85    
86 gross 1028 jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
87 gross 781 if (Finley_noError()) {
88    
89     if (grad_data_type==FINLEY_ELEMENTS) {
90     dof_offset=0;
91     s_offset=0;
92     } else if (grad_data_type==FINLEY_FACE_ELEMENTS) {
93     dof_offset=0;
94     s_offset=0;
95     } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_1) {
96     dof_offset=0;
97     s_offset=0;
98     } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2) {
99     dof_offset=numShapes;
100     s_offset=jac->ReferenceElement->Type->numShapes;
101     }
102    
103     /* check the dimensions of data */
104    
105     if (! numSamplesEqual(grad_data,jac->ReferenceElement->numQuadNodes,elements->numElements)) {
106     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples in gradient Data object");
107     } else if (! numSamplesEqual(data,1,numNodes)) {
108     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples of input Data object");
109     } else if (jac->numDim*numComps!=getDataPointSize(grad_data)) {
110     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of components in gradient data object.");
111     } else if (!isExpanded(grad_data)) {
112     Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: expanded Data object is expected for output data.");
113     } else if (! (dof_offset+jac->ReferenceElement->Type->numShapes <= numLocalNodes)) {
114     Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
115     } else if (! (s_offset+jac->ReferenceElement->Type->numShapes <= jac->ReferenceElement->Type->numNodes)) {
116     Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: offset test failed.");
117     }
118 jgs 82 }
119     /* now we can start */
120    
121 jgs 150 if (Finley_noError()) {
122 gross 853 register dim_t e,q,l,s,n;
123     register double* data_array, *grad_data_e;
124     #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
125 gross 781 {
126     if (data_type==FINLEY_NODES) {
127     if (jac->numDim==1) {
128     #define DIM 1
129     #pragma omp for schedule(static)
130     for (e=0;e<elements->numElements;e++) {
131     grad_data_e=getSampleData(grad_data,e);
132     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
133     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
134     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
135     data_array=getSampleData(data,n);
136     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
137     for (l=0;l<numComps;l++) {
138     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
139     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
140     }
141     }
142     }
143     }
144 jgs 82
145 gross 781 #undef DIM
146     } else if (jac->numDim==2) {
147     #define DIM 2
148     #pragma omp for schedule(static)
149     for (e=0;e<elements->numElements;e++) {
150     grad_data_e=getSampleData(grad_data,e);
151     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
152     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
153     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
154     data_array=getSampleData(data,n);
155     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
156     for (l=0;l<numComps;l++) {
157     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
158     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
159     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
160     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
161     }
162     }
163     }
164     }
165     #undef DIM
166     } else if (jac->numDim==3) {
167     #define DIM 3
168 gross 853 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
169 gross 781 for (e=0;e<elements->numElements;e++) {
170 gross 853 grad_data_e=getSampleData(grad_data,e);
171 gross 781 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
172     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
173     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
174     data_array=getSampleData(data,n);
175     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
176     for (l=0;l<numComps;l++) {
177 gross 853 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
178 gross 781 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
179 gross 853 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
180 gross 781 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
181 gross 853 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
182 gross 781 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
183     }
184     }
185     }
186     }
187     #undef DIM
188     }
189 gross 532
190 gross 781 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
191     if (jac->numDim==1) {
192     #define DIM 1
193     #pragma omp for schedule(static)
194     for (e=0;e<elements->numElements;e++) {
195     grad_data_e=getSampleData(grad_data,e);
196     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
197     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
198     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
199     data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
200     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
201     for (l=0;l<numComps;l++) {
202     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
203     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
204     }
205     }
206     }
207     }
208    
209     #undef DIM
210     } else if (jac->numDim==2) {
211     #define DIM 2
212     #pragma omp for schedule(static)
213     for (e=0;e<elements->numElements;e++) {
214     grad_data_e=getSampleData(grad_data,e);
215     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
216     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
217     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
218     data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
219     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
220     for (l=0;l<numComps;l++) {
221     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
222     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
223     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
224     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
225     }
226     }
227     }
228     }
229     #undef DIM
230     } else if (jac->numDim==3) {
231     #define DIM 3
232     #pragma omp for schedule(static)
233     for (e=0;e<elements->numElements;e++) {
234     grad_data_e=getSampleData(grad_data,e);
235     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
236     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
237     n=elements->Nodes[INDEX2(dof_offset+s,e,NN)];
238     data_array=getSampleData(data,nodes->degreeOfFreedom[n]);
239     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
240     for (l=0;l<numComps;l++) {
241     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
242     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
243     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
244     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
245     grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
246     jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
247     }
248     }
249     }
250     }
251     #undef DIM
252     }
253     } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
254     if (jac->numDim==1) {
255     #define DIM 1
256     #pragma omp for schedule(static)
257     for (e=0;e<elements->numElements;e++) {
258     grad_data_e=getSampleData(grad_data,e);
259     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
260     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
261     n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
262     data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
263     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
264     for (l=0;l<numComps;l++) {
265     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
266     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
267     }
268     }
269     }
270     }
271    
272     #undef DIM
273     } else if (jac->numDim==2) {
274     #define DIM 2
275     #pragma omp for schedule(static)
276     for (e=0;e<elements->numElements;e++) {
277     grad_data_e=getSampleData(grad_data,e);
278     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
279     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
280     n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
281     data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
282     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
283     for (l=0;l<numComps;l++) {
284     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
285     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
286     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
287     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
288     }
289     }
290     }
291     }
292    
293     #undef DIM
294    
295     } else if (jac->numDim==3) {
296     #define DIM 3
297     #pragma omp for schedule(static)
298     for (e=0;e<elements->numElements;e++) {
299     grad_data_e=getSampleData(grad_data,e);
300     for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
301     for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
302     n=elements->Nodes[INDEX2(elements->ReferenceElement->Type->linearNodes[dof_offset+s],e,NN)];
303     data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[n]);
304     for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
305     for (l=0;l<numComps;l++) {
306     grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
307     jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
308     grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
309     jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
310     grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
311     jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
312     }
313     }
314     }
315     }
316     #undef DIM
317     }
318     }
319     } /* end parallel region */
320 jgs 82 }
321     }
322     /*
323     * $Log$
324 jgs 150 * Revision 1.6 2005/09/15 03:44:21 jgs
325     * Merge of development branch dev-02 back to main trunk on 2005-09-15
326     *
327     * Revision 1.5.2.1 2005/09/07 06:26:17 gross
328     * the solver from finley are put into the standalone package paso now
329     *
330 jgs 123 * Revision 1.5 2005/07/08 04:07:47 jgs
331     * Merge of development branch back to main trunk on 2005-07-08
332     *
333 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
334 jgs 97 * *** empty log message ***
335 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:48 gross
336     * some changes towards 64 integers in finley
337 jgs 82 *
338 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
339     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
340 jgs 97 *
341 jgs 82 *
342 jgs 123 *
343 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26