/[escript]/branches/doubleplusgood/dudley/src/Assemble_LumpedSystem.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/dudley/src/Assemble_LumpedSystem.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26