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

Annotation of /branches/domexper/dudley/src/Assemble_PDE.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 jgs 82 /**************************************************************/
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 phornby 2078 #include "esysUtils/blocktimer.h"
45 jgs 82 #ifdef _OPENMP
46     #include <omp.h>
47     #endif
48    
49     /**************************************************************/
50    
51 jfenwick 3224 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 jgs 82
56 jfenwick 3224 bool_t reducedIntegrationOrder = FALSE;
57     char error_msg[LenErrorMsg_MAX];
58 caltinay 3247 Dudley_Assemble_Parameters p;
59 jfenwick 3224 double time0;
60     dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
61     type_t funcspace;
62     double blocktimer_start = blocktimer_time();
63 gross 798
64 jfenwick 3224 Dudley_resetError();
65 jgs 82
66 jfenwick 3224 {
67 jfenwick 3231 #ifdef ESYS_MPI
68 jfenwick 3224 int iam, numCPUs;
69 jfenwick 3231 iam = elements->MPIInfo->rank;
70     numCPUs = elements->MPIInfo->size;
71 ksteube 1312 #endif
72 jfenwick 3224 }
73 ksteube 1312
74 jfenwick 3224 if (nodes == NULL || elements == NULL)
75     return;
76     if (S == NULL && isEmpty(F))
77     return;
78 jgs 82
79 jfenwick 3224 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 jgs 82
85 jfenwick 3224 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 gross 798
90 jfenwick 3224 /* 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 jgs 82
101 jfenwick 3224 /* 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 jgs 82
129 jfenwick 3224 /* 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 jgs 82
153 jfenwick 3224 /* set all parameters in p */
154 caltinay 3247 Dudley_Assemble_getAssembleParameters(nodes, elements, S, F, reducedIntegrationOrder, &p);
155 jfenwick 3224 if (!Dudley_noError())
156     return;
157 gross 798
158 jfenwick 3224 /* check if all function spaces are the same */
159 gross 798
160 jfenwick 3224 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 jgs 82
167 jfenwick 3224 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 jgs 82
174 jfenwick 3224 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 jgs 82
181 jfenwick 3224 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 jgs 82
188 jfenwick 3224 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 jgs 82
195 jfenwick 3224 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 jgs 82
202 jfenwick 3224 /* check the dimensions: */
203 matt 1352
204 jfenwick 3224 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 jgs 82 }
259 jfenwick 3224 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 jgs 82 }
330 jfenwick 3224 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 jgs 82 }
373 jfenwick 3224 blocktimer_increment("Dudley_Assemble_PDE()", blocktimer_start);
374 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26