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

Contents of /trunk-mpi-branch/finley/src/Assemble_gradient.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1223 - (show annotations)
Fri Aug 3 02:40:39 2007 UTC (12 years, 10 months ago) by gross
File MIME type: text/plain
File size: 8855 byte(s)
first attemt towards an improved MPI version.  

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 dim_t numNodes, numShapes, numLocalNodes, numComps, NN;
37 type_t data_type=getFunctionSpaceType(data);
38 bool_t reducedShapefunction=FALSE, reducedIntegrationOrder=FALSE;
39 Finley_NodeMapping *mapping=NULL;
40 index_t dof_offset, s_offset;
41 Finley_ElementFile_Jacobeans* jac=NULL;
42 type_t grad_data_type=getFunctionSpaceType(grad_data);
43 Finley_resetError();
44 if (nodes==NULL || elements==NULL) return;
45 numComps=getDataPointSize(data);
46 NN=elements->ReferenceElement->Type->numNodes;
47 reducedIntegrationOrder=Finley_Assemble_reducedIntegrationOrder(grad_data);
48
49 if (data_type==FINLEY_NODES) {
50 reducedShapefunction=FALSE;
51 mapping=nodes->nodesMapping;
52 } else if (data_type==FINLEY_REDUCED_NODES) {
53 reducedShapefunction=TRUE;
54 mapping=nodes->reducedNodesMapping;
55 } else if (data_type==FINLEY_DEGREES_OF_FREEDOM) {
56 if (elements->MPIInfo->size>1) {
57 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: for more than one processor DEGREES_OF_FREEDOM data are not accepted as input.");
58 return;
59 }
60 reducedShapefunction=FALSE;
61 mapping=nodes->degreesOfFreedomMapping;
62 } else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) {
63 if (elements->MPIInfo->size>1) {
64 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: for more than one processor REDUCED_DEGREES_OF_FREEDOM data are not accepted as input.");
65 return;
66 }
67 reducedShapefunction=TRUE;
68 mapping=nodes->reducedDegreesOfFreedomMapping;
69 } else {
70 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.");
71 }
72
73 jac=Finley_ElementFile_borrowJacobeans(elements,nodes,reducedShapefunction,reducedIntegrationOrder);
74 if (Finley_noError()) {
75
76 numShapes=jac->ReferenceElement->Type->numShapes;
77 numLocalNodes=jac->ReferenceElement->Type->numNodes;
78 if (grad_data_type==FINLEY_CONTACT_ELEMENTS_2 || grad_data_type== FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
79 dof_offset=numShapes;
80 s_offset=jac->ReferenceElement->Type->numShapes;
81 } else {
82 dof_offset=0;
83 s_offset=0;
84 }
85
86 /* check the dimensions of data */
87
88 if (! numSamplesEqual(grad_data,jac->ReferenceElement->numQuadNodes,elements->numElements)) {
89 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples in gradient Data object");
90 } else if (! numSamplesEqual(data,1,mapping->numTargets)) {
91 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of samples of input Data object");
92 } else if (jac->numDim*numComps!=getDataPointSize(grad_data)) {
93 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: illegal number of components in gradient data object.");
94 } else if (!isExpanded(grad_data)) {
95 Finley_setError(TYPE_ERROR,"Finley_Assemble_gradient: expanded Data object is expected for output data.");
96 } else if (! (dof_offset+jac->ReferenceElement->Type->numShapes <= numLocalNodes)) {
97 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: nodes per element is inconsistent with number of jacobeans.");
98 } else if (! (s_offset+jac->ReferenceElement->Type->numShapes <= jac->ReferenceElement->Type->numNodes)) {
99 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_gradient: offset test failed.");
100 }
101 }
102 /* now we can start */
103
104 if (Finley_noError()) {
105 register dim_t e,q,l,s,n;
106 register double* data_array, *grad_data_e;
107 #pragma omp parallel private(e,q,l,s,n,data_array,grad_data_e)
108 {
109 if (jac->numDim==1) {
110 #define DIM 1
111 #pragma omp for schedule(static)
112 for (e=0;e<elements->numElements;e++) {
113 grad_data_e=getSampleData(grad_data,e);
114 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
115 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
116 n=mapping->target[elements->Nodes[INDEX2(dof_offset+s,e,NN)]];
117 data_array=getSampleData(data,n);
118 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
119 for (l=0;l<numComps;l++) {
120 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
121 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
122 }
123 }
124 }
125 }
126
127 #undef DIM
128 } else if (jac->numDim==2) {
129 #define DIM 2
130 #pragma omp for schedule(static)
131 for (e=0;e<elements->numElements;e++) {
132 grad_data_e=getSampleData(grad_data,e);
133 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
134 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
135 n=mapping->target[elements->Nodes[INDEX2(dof_offset+s,e,NN)]];
136 data_array=getSampleData(data,n);
137 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
138 for (l=0;l<numComps;l++) {
139 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
140 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
141 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
142 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
143 }
144 }
145 }
146 }
147 #undef DIM
148 } else if (jac->numDim==3) {
149 #define DIM 3
150 #pragma omp for private(e,grad_data_e,s,n,data_array,q,l) schedule(static)
151 for (e=0;e<elements->numElements;e++) {
152 grad_data_e=getSampleData(grad_data,e);
153 for (q=0;q<DIM*(jac->ReferenceElement->numQuadNodes)*numComps; q++) grad_data_e[q]=0;
154 for (s=0;s<jac->ReferenceElement->Type->numShapes;s++) {
155 n=mapping->target[elements->Nodes[INDEX2(dof_offset+s,e,NN)]];
156 data_array=getSampleData(data,n);
157 for (q=0;q<jac->ReferenceElement->numQuadNodes;q++) {
158 for (l=0;l<numComps;l++) {
159 grad_data_e[INDEX3(l,0,q,numComps,DIM)]+=data_array[l]*
160 jac->DSDX[INDEX4(s_offset+s,0,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
161 grad_data_e[INDEX3(l,1,q,numComps,DIM)]+=data_array[l]*
162 jac->DSDX[INDEX4(s_offset+s,1,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
163 grad_data_e[INDEX3(l,2,q,numComps,DIM)]+=data_array[l]*
164 jac->DSDX[INDEX4(s_offset+s,2,q,e,jac->ReferenceElement->Type->numNodes,DIM,jac->ReferenceElement->numQuadNodes)];
165 }
166 }
167 }
168 }
169 #undef DIM
170 }
171 } /* end parallel region */
172 }
173 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26