/[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 2750 - (show annotations)
Tue Nov 17 23:39:20 2009 UTC (9 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 15279 byte(s)
Fixed some compilation errors with OpenMP introduced in last commit.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 /* Disabled until the tests pass */
33 /* #define NEW_LUMPING */ /* */
34
35 /**************************************************************/
36
37 void Finley_Assemble_LumpedSystem(Finley_NodeFile* nodes,Finley_ElementFile* elements, escriptDataC* lumpedMat, escriptDataC* D)
38 {
39
40 bool_t reducedIntegrationOrder=FALSE, expandedD;
41 char error_msg[LenErrorMsg_MAX];
42 Assemble_Parameters p;
43 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s, isub;
44 type_t funcspace;
45 index_t color,*row_index=NULL;
46 __const double *D_p=NULL;
47 double *S=NULL, *EM_lumpedMat=NULL, *Vol=NULL, *lumpedMat_p=NULL;
48 register double rtmp;
49 size_t len_EM_lumpedMat_size;
50 #ifdef NEW_LUMPING
51 register double m_t=0., diagS=0.;
52 #endif
53
54 Finley_resetError();
55
56 if (nodes==NULL || elements==NULL) return;
57 if (isEmpty(lumpedMat) || isEmpty(D)) return;
58 if (isEmpty(lumpedMat) && !isEmpty(D)) {
59 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficients are non-zero but no lumped matrix is given.");
60 return;
61 }
62 funcspace=getFunctionSpaceType(D);
63 /* check if all function spaces are the same */
64 if (funcspace==FINLEY_ELEMENTS) {
65 reducedIntegrationOrder=FALSE;
66 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
67 reducedIntegrationOrder=FALSE;
68 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
69 reducedIntegrationOrder=TRUE;
70 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
71 reducedIntegrationOrder=TRUE;
72 } else {
73 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: assemblage failed because of illegal function space.");
74 }
75 if (! Finley_noError()) return;
76
77 /* set all parameters in p*/
78 Assemble_getAssembleParameters(nodes,elements,NULL,lumpedMat, reducedIntegrationOrder, &p);
79 if (! Finley_noError()) return;
80
81 /* check if all function spaces are the same */
82
83 if (! numSamplesEqual(D,p.numQuadSub,elements->numElements) ) {
84 sprintf(error_msg,"Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)",p.numQuadSub,elements->numElements);
85 Finley_setError(TYPE_ERROR,error_msg);
86 }
87
88 /* check the dimensions: */
89
90 if (p.numEqu==1) {
91 if (!isEmpty(D)) {
92 if (!isDataPointShapeEqual(D,0,dimensions)) {
93 Finley_setError(TYPE_ERROR,"Assemble_LumpedSystem: coefficient D, rank 0 expected.");
94 }
95
96 }
97 } else {
98 if (!isEmpty(D)) {
99 dimensions[0]=p.numEqu;
100 if (!isDataPointShapeEqual(D,1,dimensions)) {
101 sprintf(error_msg,"Assemble_LumpedSystem: coefficient D, expected shape (%d,)",dimensions[0]);
102 Finley_setError(TYPE_ERROR,error_msg);
103 }
104 }
105 }
106
107 if (Finley_noError()) {
108 void* buffer=allocSampleBuffer(D);
109 requireWrite(lumpedMat);
110 lumpedMat_p=getSampleDataRW(lumpedMat,0);
111 len_EM_lumpedMat=p.row_numShapesTotal*p.numEqu;
112 len_EM_lumpedMat_size=len_EM_lumpedMat*sizeof(double);
113
114 expandedD=isExpanded(D);
115 S=p.row_jac->BasisFunctions->S;
116
117 #ifdef NEW_LUMPING
118 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp, diagS, m_t)
119 #else
120 #pragma omp parallel private(color, EM_lumpedMat, row_index, Vol, D_p, s, q, k, rtmp)
121 #endif
122 {
123 EM_lumpedMat=THREAD_MEMALLOC(len_EM_lumpedMat,double);
124 row_index=THREAD_MEMALLOC(p.row_numShapesTotal,index_t);
125 if ( !Finley_checkPtr(EM_lumpedMat) && !Finley_checkPtr(row_index) ) {
126 if (p.numEqu == 1) {
127 if (expandedD) {
128
129 for (color=elements->minColor;color<=elements->maxColor;color++) {
130 /* open loop over all elements: */
131 #pragma omp for private(e) schedule(static)
132 for(e=0;e<elements->numElements;e++){
133
134 if (elements->Color[e]==color) {
135 for (isub=0; isub<p.numSub; isub++) {
136 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
137 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
138
139 D_p=getSampleDataRO(D,e,buffer);
140 #ifdef NEW_LUMPING /* HRZ lumping */
141 m_t=0; /* mass of the element: m_t */
142 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub) ];
143
144 diagS=0; /* diagonal sum: S */
145 for (s=0;s<p.row_numShapes;s++) {
146 rtmp=0;
147 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX2(q, isub,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
148 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
149 diagS+=rtmp;
150 }
151 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
152 rtmp=m_t/diagS;
153 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
154
155 #else /* row-sum lumping */
156 for (s=0;s<p.row_numShapes;s++) {
157 rtmp=0;
158 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*D_p[INDEX2(q, isub,p.numQuadSub)];
159 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp;
160 }
161 #endif
162 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)]];
163 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
164 } /* end of isub loop */
165 } /* end color check */
166
167 } /* end element loop */
168 } /* end color loop */
169 } else {
170 for (color=elements->minColor;color<=elements->maxColor;color++) {
171 /* open loop over all elements: */
172 #pragma omp for private(e) schedule(static)
173 for(e=0;e<elements->numElements;e++){
174
175 if (elements->Color[e]==color) {
176 for (isub=0; isub<p.numSub; isub++) {
177
178 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
179 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
180
181 D_p=getSampleDataRO(D,e,buffer);
182 #ifdef NEW_LUMPING /* HRZ lumping */
183 m_t=0; /* mass of the element: m_t */
184 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q];
185 m_t*=D_p[0];
186
187 diagS=0; /* diagonal sum: S */
188 for (s=0;s<p.row_numShapes;s++) {
189 rtmp=0;
190 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
191 rtmp*=D_p[0];
192 EM_lumpedMat[INDEX2(0,s,p.numEqu)]=rtmp;
193 diagS+=rtmp;
194 }
195 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
196 rtmp=m_t/diagS;
197 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(0,s,p.numEqu)]*=rtmp;
198
199 #else /* row-sum lumping */
200 for (s=0;s<p.row_numShapes;s++) {
201 rtmp=0;
202 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
203 EM_lumpedMat[INDEX2(0,s,p.numEqu)]+=rtmp*D_p[0];
204 }
205 #endif
206 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)]];
207 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
208 } /* end of isub loop */
209 } /* end color check */
210 } /* end element loop */
211 } /* end color loop */
212
213 }
214 } else {
215 if (expandedD) {
216 for (color=elements->minColor;color<=elements->maxColor;color++) {
217 /* open loop over all elements: */
218 #pragma omp for private(e) schedule(static)
219 for(e=0;e<elements->numElements;e++){
220 if (elements->Color[e]==color) {
221 for (isub=0; isub<p.numSub; isub++) {
222 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
223 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
224 D_p=getSampleDataRO(D,e,buffer);
225
226 #ifdef NEW_LUMPING /* HRZ lumping */
227 for (k=0;k<p.numEqu;k++) {
228 m_t=0; /* mass of the element: m_t */
229 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
230
231 diagS=0; /* diagonal sum: S */
232 for (s=0;s<p.row_numShapes;s++) {
233 rtmp=0;
234 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
235 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
236 diagS+=rtmp;
237 }
238 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
239 rtmp=m_t/diagS;
240 for (s=0;s<p.row_numShapes;s++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
241 }
242 #else /* row-sum lumping */
243 for (s=0;s<p.row_numShapes;s++) {
244 for (k=0;k<p.numEqu;k++) {
245 rtmp=0.;
246 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)];
247 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp;
248 }
249 }
250 #endif
251 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)]];
252 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
253 } /* end of isub loop */
254 } /* end color check */
255 } /* end element loop */
256 } /* end color loop */
257 } else {
258 for (color=elements->minColor;color<=elements->maxColor;color++) {
259 /* open loop over all elements: */
260 #pragma omp for private(e) schedule(static)
261 for(e=0;e<elements->numElements;e++){
262 if (elements->Color[e]==color) {
263 for (isub=0; isub<p.numSub; isub++) {
264 Vol=&(p.row_jac->volume[INDEX3(0,isub,e, p.numQuadSub,p.numSub)]);
265 memset(EM_lumpedMat,0,len_EM_lumpedMat_size);
266 D_p=getSampleDataRO(D,e,buffer);
267
268 #ifdef NEW_LUMPING /* HRZ lumping */
269 for (k=0;k<p.numEqu;k++) {
270 m_t=0; /* mass of the element: m_t */
271 for (q=0;q<p.numQuadSub;q++) m_t+=Vol[q]*D_p[INDEX3(k,q,isub,p.numEqu,p.numQuadSub)];
272 m_t*=D_p[k];
273 diagS=0; /* diagonal sum: S */
274 for (s=0;s<p.row_numShapes;s++) {
275 rtmp=0;
276 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)]*S[INDEX2(s,q,p.row_numShapes)];
277 rtmp*=D_p[k];
278 EM_lumpedMat[INDEX2(k,s,p.numEqu)]=rtmp;
279 diagS+=rtmp;
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++) EM_lumpedMat[INDEX2(k,s,p.numEqu)]*=rtmp;
284 }
285 #else /* row-sum lumping */
286 for (s=0;s<p.row_numShapes;s++) {
287 for (k=0;k<p.numEqu;k++) {
288 rtmp=0.;
289 for (q=0;q<p.numQuadSub;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_numShapes)];
290 EM_lumpedMat[INDEX2(k,s,p.numEqu)]+=rtmp*D_p[k];
291 }
292 }
293 #endif
294 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)]];
295 Finley_Util_AddScatter(p.row_numShapesTotal,row_index,p.numEqu,EM_lumpedMat,lumpedMat_p, p.row_DOF_UpperBound);
296 } /* end of isub loop */
297 } /* end color check */
298 } /* end element loop */
299 } /* end color loop */
300 }
301 }
302 } /* end of pointer check */
303 THREAD_MEMFREE(EM_lumpedMat);
304 THREAD_MEMFREE(row_index);
305 } /* end parallel region */
306 freeSampleBuffer(buffer);
307 }
308 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26