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

Contents of /trunk/dudley/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, 5 months ago) by caltinay
File MIME type: text/plain
File size: 11962 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 /* assembles the mass matrix in lumped form */
17
18 /* The coefficient D has to be defined on the integration points or not present. */
19
20 /* lumpedMat has to be initialized before the routine is called. */
21
22 /**************************************************************/
23
24 #include "Assemble.h"
25 #include "Util.h"
26 #ifdef _OPENMP
27 #include <omp.h>
28 #endif
29
30 #include "ShapeTable.h"
31
32 /**************************************************************/
33
34 void Dudley_Assemble_LumpedSystem(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, escriptDataC * lumpedMat,
35 escriptDataC * D, const bool_t useHRZ)
36 {
37
38 bool_t reducedIntegrationOrder = FALSE, expandedD;
39 char error_msg[LenErrorMsg_MAX];
40 Dudley_Assemble_Parameters p;
41 dim_t dimensions[ESCRIPT_MAX_DATA_RANK], k, e, len_EM_lumpedMat, q, s;
42 type_t funcspace;
43 index_t color, *row_index = NULL;
44 __const double *D_p = NULL;
45 const double *S = NULL;
46 double *EM_lumpedMat = NULL, *lumpedMat_p = NULL;
47 register double rtmp;
48 register double m_t = 0., diagS = 0.;
49
50
51 Dudley_resetError();
52
53 if (nodes == NULL || elements == NULL)
54 return;
55 if (isEmpty(lumpedMat) || isEmpty(D))
56 return;
57 if (isEmpty(lumpedMat) && !isEmpty(D))
58 {
59 Dudley_setError(TYPE_ERROR, "Dudley_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 == DUDLEY_ELEMENTS)
65 {
66 reducedIntegrationOrder = FALSE;
67 }
68 else if (funcspace == DUDLEY_FACE_ELEMENTS)
69 {
70 reducedIntegrationOrder = FALSE;
71 }
72 else if (funcspace == DUDLEY_REDUCED_ELEMENTS)
73 {
74 reducedIntegrationOrder = TRUE;
75 }
76 else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS)
77 {
78 reducedIntegrationOrder = TRUE;
79 }
80 else
81 {
82 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: assemblage failed because of illegal function space.");
83 }
84 if (!Dudley_noError())
85 return;
86
87 /* set all parameters in p */
88 Dudley_Assemble_getAssembleParameters(nodes, elements, NULL, lumpedMat, reducedIntegrationOrder, &p);
89 if (!Dudley_noError())
90 return;
91
92 /* check if all function spaces are the same */
93 if (!numSamplesEqual(D, p.numQuad, elements->numElements))
94 {
95 sprintf(error_msg, "Dudley_Assemble_LumpedSystem: sample points of coefficient D don't match (%d,%d)", p.numQuad,
96 elements->numElements);
97 Dudley_setError(TYPE_ERROR, error_msg);
98 }
99
100 /* check the dimensions: */
101 if (p.numEqu == 1)
102 {
103 if (!isEmpty(D))
104 {
105 if (!isDataPointShapeEqual(D, 0, dimensions))
106 {
107 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: coefficient D, rank 0 expected.");
108 }
109
110 }
111 }
112 else
113 {
114 if (!isEmpty(D))
115 {
116 dimensions[0] = p.numEqu;
117 if (!isDataPointShapeEqual(D, 1, dimensions))
118 {
119 sprintf(error_msg, "Dudley_Assemble_LumpedSystem: coefficient D, expected shape (%d,)", dimensions[0]);
120 Dudley_setError(TYPE_ERROR, error_msg);
121 }
122 }
123 }
124 if (Dudley_noError())
125 {
126 requireWrite(lumpedMat);
127 lumpedMat_p = getSampleDataRW(lumpedMat, 0);
128 len_EM_lumpedMat = p.numShapes * p.numEqu;
129
130 expandedD = isExpanded(D);
131 if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &S))
132 {
133 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_LumpedSystem: Unable to locate shape function.");
134 }
135 #pragma omp parallel private(color, EM_lumpedMat, row_index, D_p, s, q, k, rtmp, diagS, m_t)
136 {
137 EM_lumpedMat = THREAD_MEMALLOC(len_EM_lumpedMat, double);
138 row_index = THREAD_MEMALLOC(p.numShapes, index_t);
139 if (!Dudley_checkPtr(EM_lumpedMat) && !Dudley_checkPtr(row_index))
140 {
141 if (p.numEqu == 1)
142 { /* single equation */
143 if (expandedD)
144 { /* with expanded D */
145 for (color = elements->minColor; color <= elements->maxColor; color++)
146 {
147 /* open loop over all elements: */
148 #pragma omp for private(e) schedule(static)
149 for (e = 0; e < elements->numElements; e++)
150 {
151 if (elements->Color[e] == color)
152 {
153 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
154 D_p = getSampleDataRO(D, e);
155 if (useHRZ) {
156 m_t = 0; /* mass of the element: m_t */
157 for (q = 0; q < p.numQuad; q++)
158 m_t += vol * D_p[INDEX2(q, 0, p.numQuad)];
159 diagS = 0; /* diagonal sum: S */
160 for (s = 0; s < p.numShapes; s++)
161 {
162 rtmp = 0;
163 for (q = 0; q < p.numQuad; q++)
164 rtmp +=
165 vol * D_p[INDEX2(q, 0, p.numQuad)] * S[INDEX2(s, q, p.numShapes)] *
166 S[INDEX2(s, q, p.numShapes)];
167 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
168 diagS += rtmp;
169 }
170 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
171 rtmp = m_t / diagS;
172 for (s = 0; s < p.numShapes; s++)
173 EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
174
175 } else {/* row-sum lumping */
176 for (s = 0; s < p.numShapes; s++)
177 {
178 rtmp = 0;
179 for (q = 0; q < p.numQuad; q++)
180 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * D_p[INDEX2(q, 0, p.numQuad)];
181 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
182 }
183 }
184 for (q = 0; q < p.numShapes; q++)
185 {
186 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
187 }
188 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
189 p.row_DOF_UpperBound);
190 } /* end color check */
191 } /* end element loop */
192 } /* end color loop */
193 }
194 else
195 { /* with constant D */
196
197 for (color = elements->minColor; color <= elements->maxColor; color++)
198 {
199 /* open loop over all elements: */
200 #pragma omp for private(e) schedule(static)
201 for (e = 0; e < elements->numElements; e++)
202 {
203 if (elements->Color[e] == color)
204 {
205 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
206 D_p = getSampleDataRO(D, e);
207 if (useHRZ) { /* HRZ lumping */
208 m_t = 0; /* mass of the element: m_t */
209 for (q = 0; q < p.numQuad; q++)
210 m_t += vol;
211 diagS = 0; /* diagonal sum: S */
212 for (s = 0; s < p.numShapes; s++)
213 {
214 rtmp = 0;
215 for (q = 0; q < p.numQuad; q++)
216 {
217 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
218 }
219 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp;
220 diagS += rtmp;
221 }
222 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
223 rtmp = m_t / diagS * D_p[0];
224 for (s = 0; s < p.numShapes; s++)
225 EM_lumpedMat[INDEX2(0, s, p.numEqu)] *= rtmp;
226 } else { /* row-sum lumping */
227 for (s = 0; s < p.numShapes; s++)
228 {
229 rtmp = 0;
230 for (q = 0; q < p.numQuad; q++)
231 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
232 EM_lumpedMat[INDEX2(0, s, p.numEqu)] = rtmp * D_p[0];
233 }
234 }
235 for (q = 0; q < p.numShapes; q++)
236 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
237 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
238 p.row_DOF_UpperBound);
239 } /* end color check */
240 } /* end element loop */
241 } /* end color loop */
242
243 }
244 }
245 else
246 { /* system of equation */
247 if (expandedD)
248 { /* with expanded D */
249 for (color = elements->minColor; color <= elements->maxColor; color++)
250 {
251 /* open loop over all elements: */
252 #pragma omp for private(e) schedule(static)
253 for (e = 0; e < elements->numElements; e++)
254 {
255 if (elements->Color[e] == color)
256 {
257 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
258 D_p = getSampleDataRO(D, e);
259
260 if (useHRZ) { /* HRZ lumping */
261 for (k = 0; k < p.numEqu; k++)
262 {
263 m_t = 0; /* mass of the element: m_t */
264 for (q = 0; q < p.numQuad; q++)
265 m_t += vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
266
267 diagS = 0; /* diagonal sum: S */
268 for (s = 0; s < p.numShapes; s++)
269 {
270 rtmp = 0;
271 for (q = 0; q < p.numQuad; q++)
272 rtmp +=
273 vol * D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)] *
274 S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
275 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
276 diagS += rtmp;
277 }
278 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
279 rtmp = m_t / diagS;
280 for (s = 0; s < p.numShapes; s++)
281 EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp;
282 }
283 } else { /* row-sum lumping */
284 for (s = 0; s < p.numShapes; s++)
285 {
286 for (k = 0; k < p.numEqu; k++)
287 {
288 rtmp = 0.;
289 for (q = 0; q < p.numQuad; q++)
290 rtmp +=
291 vol * S[INDEX2(s, q, p.numShapes)] *
292 D_p[INDEX3(k, q, 0, p.numEqu, p.numQuad)];
293 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
294 }
295 }
296 }
297 for (q = 0; q < p.numShapes; q++)
298 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
299 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
300 p.row_DOF_UpperBound);
301 } /* end color check */
302 } /* end element loop */
303 } /* end color loop */
304 }
305 else
306 { /* with constant D */
307 for (color = elements->minColor; color <= elements->maxColor; color++)
308 {
309 /* open loop over all elements: */
310 #pragma omp for private(e) schedule(static)
311 for (e = 0; e < elements->numElements; e++)
312 {
313 if (elements->Color[e] == color)
314 {
315 double vol = p.row_jac->absD[e] * p.row_jac->quadweight;
316 D_p = getSampleDataRO(D, e);
317
318 if (useHRZ) { /* HRZ lumping */
319 m_t = 0; /* mass of the element: m_t */
320 for (q = 0; q < p.numQuad; q++)
321 m_t += vol;
322 diagS = 0; /* diagonal sum: S */
323 for (s = 0; s < p.numShapes; s++)
324 {
325 rtmp = 0;
326 for (q = 0; q < p.numQuad; q++)
327 rtmp += vol * S[INDEX2(s, q, p.numShapes)] * S[INDEX2(s, q, p.numShapes)];
328 for (k = 0; k < p.numEqu; k++)
329 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp;
330 diagS += rtmp;
331 }
332 /* rescale diagonals by m_t/diagS to ensure consistent mass over element */
333 rtmp = m_t / diagS;
334 for (s = 0; s < p.numShapes; s++)
335 {
336 for (k = 0; k < p.numEqu; k++)
337 EM_lumpedMat[INDEX2(k, s, p.numEqu)] *= rtmp * D_p[k];
338 }
339 } else { /* row-sum lumping */
340 for (s = 0; s < p.numShapes; s++)
341 {
342 for (k = 0; k < p.numEqu; k++)
343 {
344 rtmp = 0.;
345 for (q = 0; q < p.numQuad; q++)
346 rtmp += vol * S[INDEX2(s, q, p.numShapes)];
347 EM_lumpedMat[INDEX2(k, s, p.numEqu)] = rtmp * D_p[k];
348 }
349 }
350 }
351 for (q = 0; q < p.numShapes; q++)
352 row_index[q] = p.row_DOF[elements->Nodes[INDEX2(q, e, p.NN)]];
353 Dudley_Util_AddScatter(p.numShapes, row_index, p.numEqu, EM_lumpedMat, lumpedMat_p,
354 p.row_DOF_UpperBound);
355 } /* end color check */
356 } /* end element loop */
357 } /* end color loop */
358 }
359 }
360 } /* end of pointer check */
361 THREAD_MEMFREE(EM_lumpedMat);
362 THREAD_MEMFREE(row_index);
363 } /* end parallel region */
364 }
365 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26