← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/Assemble_PDE.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_PDE.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 system of numEq PDEs into the stiffness matrix S and right hand side F */
 
17
 
 
18
/*     -div(A*grad u)-div(B*u)+C*grad u + D*u= -div X + Y */
 
19
 
 
20
/*      -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m = -(X_{k,i})_i + Y_k */
 
21
 
 
22
/*    u has numComp components. */
 
23
 
 
24
/*    Shape of the coefficients: */
 
25
 
 
26
/*      A = numEqu x numDim x numComp x numDim */
 
27
/*      B = numDim x numEqu x numComp  */
 
28
/*      C = numEqu x numDim x numComp  */
 
29
/*      D = numEqu x numComp  */
 
30
/*      X = numEqu x numDim   */
 
31
/*      Y = numEqu */
 
32
 
 
33
/*    The coefficients A,B,C,D,X and Y have to be defined on the integartion points or not present (=NULL). */
 
34
 
 
35
/*    S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
 
36
/*    the right hand side of the PDE is not processed.  */
 
37
 
 
38
/*    The routine does not consider any boundary conditions. */
 
39
 
 
40
/**************************************************************/
 
41
 
 
42
#include "Assemble.h"
 
43
#include "Util.h"
 
44
#include "esysUtils/blocktimer.h"
 
45
#ifdef _OPENMP
 
46
#include <omp.h>
 
47
#endif
 
48
 
 
49
/**************************************************************/
 
50
 
 
51
void Dudley_Assemble_PDE(Dudley_NodeFile * nodes, Dudley_ElementFile * elements, Paso_SystemMatrix * S,
 
52
                         escriptDataC * F, escriptDataC * A, escriptDataC * B, escriptDataC * C, escriptDataC * D,
 
53
                         escriptDataC * X, escriptDataC * Y)
 
54
{
 
55
 
 
56
    bool_t reducedIntegrationOrder = FALSE;
 
57
    char error_msg[LenErrorMsg_MAX];
 
58
    Dudley_Assemble_Parameters p;
 
59
    double time0;
 
60
    dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
 
61
    type_t funcspace;
 
62
    double blocktimer_start = blocktimer_time();
 
63
 
 
64
    Dudley_resetError();
 
65
 
 
66
    {
 
67
#ifdef ESYS_MPI
 
68
        int iam, numCPUs;
 
69
        iam = elements->MPIInfo->rank;
 
70
        numCPUs = elements->MPIInfo->size;
 
71
#endif
 
72
    }
 
73
 
 
74
    if (nodes == NULL || elements == NULL)
 
75
        return;
 
76
    if (S == NULL && isEmpty(F))
 
77
        return;
 
78
 
 
79
    if (isEmpty(F) && !isEmpty(X) && !isEmpty(F))
 
80
    {
 
81
        Dudley_setError(TYPE_ERROR,
 
82
                        "Dudley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
 
83
    }
 
84
 
 
85
    if (S == NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D))
 
86
    {
 
87
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
 
88
    }
 
89
 
 
90
    /*  get the functionspace for this assemblage call */
 
91
    funcspace = UNKNOWN;
 
92
    updateFunctionSpaceType(funcspace, A);
 
93
    updateFunctionSpaceType(funcspace, B);
 
94
    updateFunctionSpaceType(funcspace, C);
 
95
    updateFunctionSpaceType(funcspace, D);
 
96
    updateFunctionSpaceType(funcspace, X);
 
97
    updateFunctionSpaceType(funcspace, Y);
 
98
    if (funcspace == UNKNOWN)
 
99
        return;                 /* all  data are empty */
 
100
 
 
101
    /* check if all function spaces are the same */
 
102
    if (!functionSpaceTypeEqual(funcspace, A))
 
103
    {
 
104
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient A");
 
105
    }
 
106
    if (!functionSpaceTypeEqual(funcspace, B))
 
107
    {
 
108
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient B");
 
109
    }
 
110
    if (!functionSpaceTypeEqual(funcspace, C))
 
111
    {
 
112
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient C");
 
113
    }
 
114
    if (!functionSpaceTypeEqual(funcspace, D))
 
115
    {
 
116
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient D");
 
117
    }
 
118
    if (!functionSpaceTypeEqual(funcspace, X))
 
119
    {
 
120
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient X");
 
121
    }
 
122
    if (!functionSpaceTypeEqual(funcspace, Y))
 
123
    {
 
124
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient Y");
 
125
    }
 
126
    if (!Dudley_noError())
 
127
        return;
 
128
 
 
129
    /* check if all function spaces are the same */
 
130
    if (funcspace == DUDLEY_ELEMENTS)
 
131
    {
 
132
        reducedIntegrationOrder = FALSE;
 
133
    }
 
134
    else if (funcspace == DUDLEY_FACE_ELEMENTS)
 
135
    {
 
136
        reducedIntegrationOrder = FALSE;
 
137
    }
 
138
    else if (funcspace == DUDLEY_REDUCED_ELEMENTS)
 
139
    {
 
140
        reducedIntegrationOrder = TRUE;
 
141
    }
 
142
    else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS)
 
143
    {
 
144
        reducedIntegrationOrder = TRUE;
 
145
    }
 
146
    else
 
147
    {
 
148
        Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: assemblage failed because of illegal function space.");
 
149
    }
 
150
    if (!Dudley_noError())
 
151
        return;
 
152
 
 
153
    /* set all parameters in p */
 
154
    Dudley_Assemble_getAssembleParameters(nodes, elements, S, F, reducedIntegrationOrder, &p);
 
155
    if (!Dudley_noError())
 
156
        return;
 
157
 
 
158
    /* check if all function spaces are the same */
 
159
 
 
160
    if (!numSamplesEqual(A, p.numQuad, elements->numElements))
 
161
    {
 
162
        sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)", p.numQuad,
 
163
                elements->numElements);
 
164
        Dudley_setError(TYPE_ERROR, error_msg);
 
165
    }
 
166
 
 
167
    if (!numSamplesEqual(B, p.numQuad, elements->numElements))
 
168
    {
 
169
        sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)", p.numQuad,
 
170
                elements->numElements);
 
171
        Dudley_setError(TYPE_ERROR, error_msg);
 
172
    }
 
173
 
 
174
    if (!numSamplesEqual(C, p.numQuad, elements->numElements))
 
175
    {
 
176
        sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)", p.numQuad,
 
177
                elements->numElements);
 
178
        Dudley_setError(TYPE_ERROR, error_msg);
 
179
    }
 
180
 
 
181
    if (!numSamplesEqual(D, p.numQuad, elements->numElements))
 
182
    {
 
183
        sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)", p.numQuad,
 
184
                elements->numElements);
 
185
        Dudley_setError(TYPE_ERROR, error_msg);
 
186
    }
 
187
 
 
188
    if (!numSamplesEqual(X, p.numQuad, elements->numElements))
 
189
    {
 
190
        sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)", p.numQuad,
 
191
                elements->numElements);
 
192
        Dudley_setError(TYPE_ERROR, error_msg);
 
193
    }
 
194
 
 
195
    if (!numSamplesEqual(Y, p.numQuad, elements->numElements))
 
196
    {
 
197
        sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)", p.numQuad,
 
198
                elements->numElements);
 
199
        Dudley_setError(TYPE_ERROR, error_msg);
 
200
    }
 
201
 
 
202
    /*  check the dimensions: */
 
203
 
 
204
    if (p.numEqu == 1 && p.numComp == 1)
 
205
    {
 
206
        if (!isEmpty(A))
 
207
        {
 
208
            dimensions[0] = p.numDim;
 
209
            dimensions[1] = p.numDim;
 
210
            if (!isDataPointShapeEqual(A, 2, dimensions))
 
211
            {
 
212
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",
 
213
                        dimensions[0], dimensions[1]);
 
214
                Dudley_setError(TYPE_ERROR, error_msg);
 
215
            }
 
216
        }
 
217
        if (!isEmpty(B))
 
218
        {
 
219
            dimensions[0] = p.numDim;
 
220
            if (!isDataPointShapeEqual(B, 1, dimensions))
 
221
            {
 
222
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient B: illegal shape (%d,)", dimensions[0]);
 
223
                Dudley_setError(TYPE_ERROR, error_msg);
 
224
            }
 
225
        }
 
226
        if (!isEmpty(C))
 
227
        {
 
228
            dimensions[0] = p.numDim;
 
229
            if (!isDataPointShapeEqual(C, 1, dimensions))
 
230
            {
 
231
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient C, expected shape (%d,)", dimensions[0]);
 
232
                Dudley_setError(TYPE_ERROR, error_msg);
 
233
            }
 
234
        }
 
235
        if (!isEmpty(D))
 
236
        {
 
237
            if (!isDataPointShapeEqual(D, 0, dimensions))
 
238
            {
 
239
                Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficient D, rank 0 expected.");
 
240
            }
 
241
        }
 
242
        if (!isEmpty(X))
 
243
        {
 
244
            dimensions[0] = p.numDim;
 
245
            if (!isDataPointShapeEqual(X, 1, dimensions))
 
246
            {
 
247
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient X, expected shape (%d,", dimensions[0]);
 
248
                Dudley_setError(TYPE_ERROR, error_msg);
 
249
            }
 
250
        }
 
251
        if (!isEmpty(Y))
 
252
        {
 
253
            if (!isDataPointShapeEqual(Y, 0, dimensions))
 
254
            {
 
255
                Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficient Y, rank 0 expected.");
 
256
            }
 
257
        }
 
258
    }
 
259
    else
 
260
    {
 
261
        if (!isEmpty(A))
 
262
        {
 
263
            dimensions[0] = p.numEqu;
 
264
            dimensions[1] = p.numDim;
 
265
            dimensions[2] = p.numComp;
 
266
            dimensions[3] = p.numDim;
 
267
            if (!isDataPointShapeEqual(A, 4, dimensions))
 
268
            {
 
269
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)", dimensions[0],
 
270
                        dimensions[1], dimensions[2], dimensions[3]);
 
271
                Dudley_setError(TYPE_ERROR, error_msg);
 
272
            }
 
273
        }
 
274
        if (!isEmpty(B))
 
275
        {
 
276
            dimensions[0] = p.numEqu;
 
277
            dimensions[1] = p.numDim;
 
278
            dimensions[2] = p.numComp;
 
279
            if (!isDataPointShapeEqual(B, 3, dimensions))
 
280
            {
 
281
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)", dimensions[0],
 
282
                        dimensions[1], dimensions[2]);
 
283
                Dudley_setError(TYPE_ERROR, error_msg);
 
284
            }
 
285
        }
 
286
        if (!isEmpty(C))
 
287
        {
 
288
            dimensions[0] = p.numEqu;
 
289
            dimensions[1] = p.numComp;
 
290
            dimensions[2] = p.numDim;
 
291
            if (!isDataPointShapeEqual(C, 3, dimensions))
 
292
            {
 
293
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)", dimensions[0],
 
294
                        dimensions[1], dimensions[2]);
 
295
                Dudley_setError(TYPE_ERROR, error_msg);
 
296
            }
 
297
        }
 
298
        if (!isEmpty(D))
 
299
        {
 
300
            dimensions[0] = p.numEqu;
 
301
            dimensions[1] = p.numComp;
 
302
            if (!isDataPointShapeEqual(D, 2, dimensions))
 
303
            {
 
304
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient D, expected shape (%d,%d)", dimensions[0],
 
305
                        dimensions[1]);
 
306
                Dudley_setError(TYPE_ERROR, error_msg);
 
307
            }
 
308
        }
 
309
        if (!isEmpty(X))
 
310
        {
 
311
            dimensions[0] = p.numEqu;
 
312
            dimensions[1] = p.numDim;
 
313
            if (!isDataPointShapeEqual(X, 2, dimensions))
 
314
            {
 
315
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient X, expected shape (%d,%d)", dimensions[0],
 
316
                        dimensions[1]);
 
317
                Dudley_setError(TYPE_ERROR, error_msg);
 
318
            }
 
319
        }
 
320
        if (!isEmpty(Y))
 
321
        {
 
322
            dimensions[0] = p.numEqu;
 
323
            if (!isDataPointShapeEqual(Y, 1, dimensions))
 
324
            {
 
325
                sprintf(error_msg, "Dudley_Assemble_PDE: coefficient Y, expected shape (%d,)", dimensions[0]);
 
326
                Dudley_setError(TYPE_ERROR, error_msg);
 
327
            }
 
328
        }
 
329
    }
 
330
    if (Dudley_noError())
 
331
    {
 
332
        time0 = Dudley_timer();
 
333
        if (p.numEqu == p.numComp)
 
334
        {
 
335
            if (p.numEqu > 1)
 
336
            {
 
337
                /* system of PDESs */
 
338
                if (p.numDim == 3)
 
339
                {
 
340
                    Dudley_Assemble_PDE_System2_3D(p, elements, S, F, A, B, C, D, X, Y);
 
341
                }
 
342
                else if (p.numDim == 2)
 
343
                {
 
344
                    Dudley_Assemble_PDE_System2_2D(p, elements, S, F, A, B, C, D, X, Y);
 
345
                }
 
346
                else
 
347
                {
 
348
                    Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
 
349
                }
 
350
            }
 
351
            else
 
352
            {
 
353
                /* single PDES */
 
354
                if (p.numDim == 3)
 
355
                {
 
356
                    Dudley_Assemble_PDE_Single2_3D(p, elements, S, F, A, B, C, D, X, Y);
 
357
                }
 
358
                else if (p.numDim == 2)
 
359
                {
 
360
                    Dudley_Assemble_PDE_Single2_2D(p, elements, S, F, A, B, C, D, X, Y);
 
361
                }
 
362
                else
 
363
                {
 
364
                    Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
 
365
                }
 
366
            }
 
367
        }
 
368
        else
 
369
        {
 
370
            Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE requires number of equations == number of solutions  .");
 
371
        }
 
372
    }
 
373
    blocktimer_increment("Dudley_Assemble_PDE()", blocktimer_start);
 
374
}