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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3490 - (show annotations)
Wed Mar 30 02:24:33 2011 UTC (8 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 14974 byte(s)
More gcc-4.6 fixes (mostly initialized-but-unused-var warnings)

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* assembles the mass matrix in lumped form */
18
19 /* The coefficient D has to be defined on the integration points or not present. */
20
21 /* lumpedMat has to be initialized before the routine is called. */
22
23 /**************************************************************/
24
25 #include "Assemble.h"
26 #include "Util.h"
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
30
31
32 /**************************************************************/
33
34 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D, const bool_t useHRZ)
35 {
36
37 bool_t reducedIntegrationOrder=FALSE, expandedD;
38 char error_msg[LenErrorMsg_MAX];
39 Finley_Assemble_Parameters p;
40 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s, isub;
41 type_t funcspace;
42 index_t color,*row_index=NULL;
43 __const double *D_p=NULL;
44 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
45 register double rtmp;
46 register double m_t=0., diagS=0., rtmp2=0.;
47
48 Finley_resetError();
49
50 if (nodes==NULL || elements==NULL) return;
51 if (isEmpty(lumpedMat) || isEmpty(D)) return;
52 if (isEmpty(lumpedMat) && !isEmpty(D)) {
53 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
54 return;
55 }
56 funcspace=getFunctionSpaceType(D);
57 /* check if all function spaces are the same */
58 if (funcspace==FINLEY_ELEMENTS) {
59 reducedIntegrationOrder=FALSE;
60 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
61 reducedIntegrationOrder=FALSE;
62 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
63 reducedIntegrationOrder=TRUE;
64 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
65 reducedIntegrationOrder=TRUE;
66 } else {
67 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: assemblage failed because of illegal function space.");
68 }
69 if (! Finley_noError()) return;
70
71 /* set all parameters in p*/
72 Finley_Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
73 if (! Finley_noError()) return;
74
75 /* check if all function spaces are the same */
76 if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
77 sprintf(error_msg,"Finley_Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
78 Finley_setError(TYPE_ERROR,error_msg);
79 }
80
81 /* check the dimensions: */
82 if (p.numEqu==1) {
83 if (!isEmpty(D)) {
84 if (!isDataPointShapeEqual(D,0,dimensions)) {
85 Finley_setError(TYPE_ERROR,"Finley_Assemble_LumpedSystem: coefficient D, rank 0 expected.");
86 }
87
88 }
89 } else {
90 if (!isEmpty(D)) {
91 dimensions[0]=p.numEqu;
92 if (!isDataPointShapeEqual(D,1,dimensions)) {
93 sprintf(error_msg,"Finley_Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
94 Finley_setError(TYPE_ERROR,error_msg);
95 }
96 }
97 }
98 if (Finley_noError()) {
99 requireWrite(lumpedMat);
100 lumpedMat_p=getSampleDataRW(lumpedMat,0);
101 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
102
103 expandedD=isExpanded(D);
104 S=p.row_jac->BasisFunctions->S;
105
106 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t, isub, rtmp2)
107 {
108 EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
109 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
110 if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
111 if (p.numEqu == 1) { /* single equation */
112 if (expandedD) { /* with expanded D */
113 for (color=elements->minColor;color<=elements->maxColor;color++) {
114 /* open loop over all elements: */
115 #pragma omp for private(e) schedule(static)
116 for(e=0;e<elements->numElements;e++){
117 if (elements->Color[e]==color) {
118 for (isub=0; isub<p.numSub; isub++) {
119 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
120 D_p=getSampleDataRO(D,e);
121 if (useHRZ) {
122
123 m_t=0; /* mass of the element: m_t */
124 #pragma ivdep
125 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];
126
127 diagS=0; /* diagonal sum: S */
128 for (s=0;s<p.row_numShapes;s++) {
129 rtmp=0;
130 #pragma ivdep
131 for (q=0;q<p.numQuadSub;q++) {
132 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
133 rtmp+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub)]*rtmp2*rtmp2;
134 }
135 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
136 diagS+=rtmp;
137 }
138 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
139 rtmp=m_t/diagS;
140 #pragma ivdep
141 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
142
143 } else { /* row-sum lumping */
144 for (s=0;s<p.row_numShapes;s++) {
145 rtmp=0;
146 #pragma ivdep
147 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];
148 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
149 }
150
151 }
152 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
153 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
154 } /* end of isub loop */
155 } /* end color check */
156 } /* end element loop */
157 } /* end color loop */
158
159
160 } else { /* with constant D */
161
162 for (color=elements->minColor;color<=elements->maxColor;color++) {
163 /* open loop over all elements: */
164 #pragma omp for private(e) schedule(static)
165 for(e=0;e<elements->numElements;e++){
166 if (elements->Color[e]==color) {
167 for (isub=0; isub<p.numSub; isub++) {
168 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
169 D_p=getSampleDataRO(D,e);
170 if (useHRZ) { /* HRZ lumping */
171 m_t=0; /* mass of the element: m_t */
172 #pragma ivdep
173 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
174 diagS=0; /* diagonal sum: S */
175 for (s=0;s<p.row_numShapes;s++) {
176 rtmp=0;
177 #pragma ivdep
178 for (q=0;q<p.numQuadSub;q++){
179 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
180 rtmp+=Vol[q]*rtmp2*rtmp2;
181 }
182 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
183 diagS+=rtmp;
184 }
185 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
186 rtmp=m_t/diagS*D_p[0];
187 #pragma ivdep
188 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
189 } else { /* row-sum lumping */
190 for (s=0;s<p.row_numShapes;s++) {
191 rtmp=0;
192 #pragma ivdep
193 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
194 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp*D_p[0];
195 }
196 }
197 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
198 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
199 } /* end of isub loop */
200 } /* end color check */
201 } /* end element loop */
202 } /* end color loop */
203
204 }
205 } else { /* system of equation */
206 if (expandedD) { /* with expanded D */
207 for (color=elements->minColor;color<=elements->maxColor;color++) {
208 /* open loop over all elements: */
209 #pragma omp for private(e) schedule(static)
210 for(e=0;e<elements->numElements;e++){
211 if (elements->Color[e]==color) {
212 for (isub=0; isub<p.numSub; isub++) {
213 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
214 D_p=getSampleDataRO(D,e);
215
216 if (useHRZ) { /* HRZ lumping */
217 for (k=0;k<p.numEqu;k++) {
218 m_t=0; /* mass of the element: m_t */
219 #pragma ivdep
220 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
221
222 diagS=0; /* diagonal sum: S */
223 for (s=0;s<p.row_numShapes;s++) {
224 rtmp=0;
225 #pragma ivdep
226 for (q=0;q<p.numQuadSub;q++) {
227 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
228 rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*rtmp2*rtmp2;
229 }
230 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
231 diagS+=rtmp;
232 }
233 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
234 rtmp=m_t/diagS;
235 #pragma ivdep
236 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
237 }
238 } else { /* row-sum lumping */
239 for (s=0;s<p.row_numShapes;s++) {
240 for (k=0;k<p.numEqu;k++) {
241 rtmp=0.;
242 #pragma ivdep
243 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
244 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
245 }
246 }
247 }
248 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
249 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
250 } /* end of isub loop */
251 } /* end color check */
252 } /* end element loop */
253 } /* end color loop */
254 } else { /* with constant D */
255 for (color=elements->minColor;color<=elements->maxColor;color++) {
256 /* open loop over all elements: */
257 #pragma omp for private(e) schedule(static)
258 for(e=0;e<elements->numElements;e++){
259 if (elements->Color[e]==color) {
260 for (isub=0; isub<p.numSub; isub++) {
261 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
262 D_p=getSampleDataRO(D,e);
263
264 if (useHRZ) { /* HRZ lumping */
265 m_t=0; /* mass of the element: m_t */
266 #pragma ivdep
267 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
268 diagS=0; /* diagonal sum: S */
269 for (s=0;s<p.row_numShapes;s++) {
270 rtmp=0;
271 #pragma ivdep
272 for (q=0;q<p.numQuadSub;q++) {
273 rtmp2=S[INDEX2(s,q,p.row_numShapes)];
274 rtmp+=Vol[q]*rtmp2*rtmp2;
275 }
276 #pragma ivdep
277 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
278 diagS+=rtmp;
279 }
280
281 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
282 rtmp=m_t/diagS;
283 for (s=0;s<p.row_numShapes;s++) {
284 #pragma ivdep
285 for (k=0;k<p.numEqu;k++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp*D_p[k];
286 }
287 } else { /* row-sum lumping */
288 for (s=0;s<p.row_numShapes;s++) {
289 for (k=0;k<p.numEqu;k++) {
290 rtmp=0.;
291 #pragma ivdep
292 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
293 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp*D_p[k];
294 }
295 }
296 }
297 for (q=0;q<p.row_numShapesTotal;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[INDEX2(q,isub,p.row_numShapesTotal)],e,p.NN)]];
298 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
299 } /* end of isub loop */
300 } /* end color check */
301 } /* end element loop */
302 } /* end color loop */
303 }
304 }
305 } /* end of pointer check */
306 THREAD_MEMFREE(EM_lumpedMat);
307 THREAD_MEMFREE(row_index);
308 } /* end parallel region */
309 }
310 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26