← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/Assemble_LumpedSystem.c

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/src/Assemble_LumpedSystem.c
 
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
}