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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26