/[escript]/branches/domexper/dudley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (8 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 11943 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26