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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 855 - (show annotations)
Fri Sep 22 09:30:06 2006 UTC (13 years, 1 month ago) by gross
File MIME type: text/plain
File size: 18109 byte(s)
some tests for triangle meshes added
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 /* assemblage rjacines: calculate the gradient of nodal data at quadrature points */
16
17 /**************************************************************/
18
19 /* Copyrights by ACcESS Australia, 2003,2004,2005 */
20 /* author: gross@access.edu.au */
21 /* version: $Id$ */
22
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 escriptDataC* grad_data,escriptDataC* data) {
35
36 Finley_resetError();
37 if (nodes==NULL || elements==NULL) return;
38 dim_t numComps=getDataPointSize(data);
39 dim_t NN=elements->ReferenceElement->Type->numNodes;
40 dim_t numNodes, numShapes, numLocalNodes;
41 index_t dof_offset, s_offset;
42 type_t data_type=getFunctionSpaceType(data);
43 type_t grad_data_type=getFunctionSpaceType(grad_data);
44 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
45
46 if (data_type==FINLEY_NODES) {
47 reducedShapefunction=FALSE;
48 numNodes=nodes->numNodes;
49 numShapes=elements->ReferenceElement->Type->numShapes;
50 numLocalNodes=elements->ReferenceElement->Type->numNodes;
51 }
52 /* lock these two options jac for the MPI version */
53 #ifndef PASO_MPI
54 else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
55 reducedShapefunction=FALSE;
56 numNodes=nodes->numDegreesOfFreedom;
57 numShapes=elements->ReferenceElement->Type->numShapes;
58 numLocalNodes=elements->ReferenceElement->Type->numNodes;
59 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
60 reducedShapefunction=TRUE;
61 numNodes=nodes->reducedNumDegreesOfFreedom;
62 numShapes=elements->LinearReferenceElement->Type->numShapes;
63 numLocalNodes=elements->LinearReferenceElement->Type->numNodes;
64 }
65 #endif
66 else {
67 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
68 }
69 if (grad_data_type==FINLEY_ELEMENTS) {
70 reducedIntegrationOrder=FALSE;
71 dof_offset=0;
72 } else if (grad_data_type==FINLEY_FACE_ELEMENTS) {
73 reducedIntegrationOrder=FALSE;
74 dof_offset=0;
75 } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_1) {
76 reducedIntegrationOrder=FALSE;
77 dof_offset=0;
78 } else if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2) {
79 reducedIntegrationOrder=FALSE;
80 /* don't reset offset */
81 } else {
82 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculated for requested location.");
83 }
84
85 Finley_ElementFile_Jacobeans* jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
86
87 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 }
119 /* now we can start */
120
121 if (Finley_noError()) {
122 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 {
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
145 #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 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
169 for (e=0;e<elements->numElements;e++) {
170 grad_data_e=getSampleData(grad_data,e);
171 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 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
178 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
179 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
180 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
181 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
182 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
190 } 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 }
321 }
322 /*
323 * $Log$
324 * 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 * 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 * Revision 1.4 2004/12/15 07:08:32 jgs
334 * *** empty log message ***
335 * Revision 1.1.1.1.2.2 2005/06/29 02:34:48 gross
336 * some changes towards 64 integers in finley
337 *
338 * 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 *
341 *
342 *
343 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26