/[escript]/trunk/dudley/src/Assemble_PDE.c
ViewVC logotype

Contents of /trunk/dudley/src/Assemble_PDE.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3490 - (show annotations)
Wed Mar 30 02:24:33 2011 UTC (8 years, 4 months ago) by caltinay
File MIME type: text/plain
File size: 10987 byte(s)
More gcc-4.6 fixes (mostly initialized-but-unused-var warnings)

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 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
60 type_t funcspace;
61 double blocktimer_start = blocktimer_time();
62
63 Dudley_resetError();
64
65 if (nodes == NULL || elements == NULL)
66 return;
67 if (S == NULL && isEmpty(F))
68 return;
69
70 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F))
71 {
72 Dudley_setError(TYPE_ERROR,
73 "Dudley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
74 }
75
76 if (S == NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D))
77 {
78 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
79 }
80
81 /* get the functionspace for this assemblage call */
82 funcspace = UNKNOWN;
83 updateFunctionSpaceType(funcspace, A);
84 updateFunctionSpaceType(funcspace, B);
85 updateFunctionSpaceType(funcspace, C);
86 updateFunctionSpaceType(funcspace, D);
87 updateFunctionSpaceType(funcspace, X);
88 updateFunctionSpaceType(funcspace, Y);
89 if (funcspace == UNKNOWN)
90 return; /* all data are empty */
91
92 /* check if all function spaces are the same */
93 if (!functionSpaceTypeEqual(funcspace, A))
94 {
95 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient A");
96 }
97 if (!functionSpaceTypeEqual(funcspace, B))
98 {
99 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient B");
100 }
101 if (!functionSpaceTypeEqual(funcspace, C))
102 {
103 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient C");
104 }
105 if (!functionSpaceTypeEqual(funcspace, D))
106 {
107 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient D");
108 }
109 if (!functionSpaceTypeEqual(funcspace, X))
110 {
111 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient X");
112 }
113 if (!functionSpaceTypeEqual(funcspace, Y))
114 {
115 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: unexpected function space type for coefficient Y");
116 }
117 if (!Dudley_noError())
118 return;
119
120 /* check if all function spaces are the same */
121 if (funcspace == DUDLEY_ELEMENTS)
122 {
123 reducedIntegrationOrder = FALSE;
124 }
125 else if (funcspace == DUDLEY_FACE_ELEMENTS)
126 {
127 reducedIntegrationOrder = FALSE;
128 }
129 else if (funcspace == DUDLEY_REDUCED_ELEMENTS)
130 {
131 reducedIntegrationOrder = TRUE;
132 }
133 else if (funcspace == DUDLEY_REDUCED_FACE_ELEMENTS)
134 {
135 reducedIntegrationOrder = TRUE;
136 }
137 else
138 {
139 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: assemblage failed because of illegal function space.");
140 }
141 if (!Dudley_noError())
142 return;
143
144 /* set all parameters in p */
145 Dudley_Assemble_getAssembleParameters(nodes, elements, S, F, reducedIntegrationOrder, &p);
146 if (!Dudley_noError())
147 return;
148
149 /* check if all function spaces are the same */
150
151 if (!numSamplesEqual(A, p.numQuad, elements->numElements))
152 {
153 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)", p.numQuad,
154 elements->numElements);
155 Dudley_setError(TYPE_ERROR, error_msg);
156 }
157
158 if (!numSamplesEqual(B, p.numQuad, elements->numElements))
159 {
160 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)", p.numQuad,
161 elements->numElements);
162 Dudley_setError(TYPE_ERROR, error_msg);
163 }
164
165 if (!numSamplesEqual(C, p.numQuad, elements->numElements))
166 {
167 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)", p.numQuad,
168 elements->numElements);
169 Dudley_setError(TYPE_ERROR, error_msg);
170 }
171
172 if (!numSamplesEqual(D, p.numQuad, elements->numElements))
173 {
174 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)", p.numQuad,
175 elements->numElements);
176 Dudley_setError(TYPE_ERROR, error_msg);
177 }
178
179 if (!numSamplesEqual(X, p.numQuad, elements->numElements))
180 {
181 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)", p.numQuad,
182 elements->numElements);
183 Dudley_setError(TYPE_ERROR, error_msg);
184 }
185
186 if (!numSamplesEqual(Y, p.numQuad, elements->numElements))
187 {
188 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)", p.numQuad,
189 elements->numElements);
190 Dudley_setError(TYPE_ERROR, error_msg);
191 }
192
193 /* check the dimensions: */
194
195 if (p.numEqu == 1 && p.numComp == 1)
196 {
197 if (!isEmpty(A))
198 {
199 dimensions[0] = p.numDim;
200 dimensions[1] = p.numDim;
201 if (!isDataPointShapeEqual(A, 2, dimensions))
202 {
203 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",
204 dimensions[0], dimensions[1]);
205 Dudley_setError(TYPE_ERROR, error_msg);
206 }
207 }
208 if (!isEmpty(B))
209 {
210 dimensions[0] = p.numDim;
211 if (!isDataPointShapeEqual(B, 1, dimensions))
212 {
213 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient B: illegal shape (%d,)", dimensions[0]);
214 Dudley_setError(TYPE_ERROR, error_msg);
215 }
216 }
217 if (!isEmpty(C))
218 {
219 dimensions[0] = p.numDim;
220 if (!isDataPointShapeEqual(C, 1, dimensions))
221 {
222 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient C, expected shape (%d,)", dimensions[0]);
223 Dudley_setError(TYPE_ERROR, error_msg);
224 }
225 }
226 if (!isEmpty(D))
227 {
228 if (!isDataPointShapeEqual(D, 0, dimensions))
229 {
230 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficient D, rank 0 expected.");
231 }
232 }
233 if (!isEmpty(X))
234 {
235 dimensions[0] = p.numDim;
236 if (!isDataPointShapeEqual(X, 1, dimensions))
237 {
238 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient X, expected shape (%d,", dimensions[0]);
239 Dudley_setError(TYPE_ERROR, error_msg);
240 }
241 }
242 if (!isEmpty(Y))
243 {
244 if (!isDataPointShapeEqual(Y, 0, dimensions))
245 {
246 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficient Y, rank 0 expected.");
247 }
248 }
249 }
250 else
251 {
252 if (!isEmpty(A))
253 {
254 dimensions[0] = p.numEqu;
255 dimensions[1] = p.numDim;
256 dimensions[2] = p.numComp;
257 dimensions[3] = p.numDim;
258 if (!isDataPointShapeEqual(A, 4, dimensions))
259 {
260 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)", dimensions[0],
261 dimensions[1], dimensions[2], dimensions[3]);
262 Dudley_setError(TYPE_ERROR, error_msg);
263 }
264 }
265 if (!isEmpty(B))
266 {
267 dimensions[0] = p.numEqu;
268 dimensions[1] = p.numDim;
269 dimensions[2] = p.numComp;
270 if (!isDataPointShapeEqual(B, 3, dimensions))
271 {
272 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)", dimensions[0],
273 dimensions[1], dimensions[2]);
274 Dudley_setError(TYPE_ERROR, error_msg);
275 }
276 }
277 if (!isEmpty(C))
278 {
279 dimensions[0] = p.numEqu;
280 dimensions[1] = p.numComp;
281 dimensions[2] = p.numDim;
282 if (!isDataPointShapeEqual(C, 3, dimensions))
283 {
284 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)", dimensions[0],
285 dimensions[1], dimensions[2]);
286 Dudley_setError(TYPE_ERROR, error_msg);
287 }
288 }
289 if (!isEmpty(D))
290 {
291 dimensions[0] = p.numEqu;
292 dimensions[1] = p.numComp;
293 if (!isDataPointShapeEqual(D, 2, dimensions))
294 {
295 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient D, expected shape (%d,%d)", dimensions[0],
296 dimensions[1]);
297 Dudley_setError(TYPE_ERROR, error_msg);
298 }
299 }
300 if (!isEmpty(X))
301 {
302 dimensions[0] = p.numEqu;
303 dimensions[1] = p.numDim;
304 if (!isDataPointShapeEqual(X, 2, dimensions))
305 {
306 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient X, expected shape (%d,%d)", dimensions[0],
307 dimensions[1]);
308 Dudley_setError(TYPE_ERROR, error_msg);
309 }
310 }
311 if (!isEmpty(Y))
312 {
313 dimensions[0] = p.numEqu;
314 if (!isDataPointShapeEqual(Y, 1, dimensions))
315 {
316 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient Y, expected shape (%d,)", dimensions[0]);
317 Dudley_setError(TYPE_ERROR, error_msg);
318 }
319 }
320 }
321 if (Dudley_noError())
322 {
323 if (p.numEqu == p.numComp)
324 {
325 if (p.numEqu > 1)
326 {
327 /* system of PDESs */
328 if (p.numDim == 3)
329 {
330 Dudley_Assemble_PDE_System2_3D(p, elements, S, F, A, B, C, D, X, Y);
331 }
332 else if (p.numDim == 2)
333 {
334 Dudley_Assemble_PDE_System2_2D(p, elements, S, F, A, B, C, D, X, Y);
335 }
336 else
337 {
338 Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
339 }
340 }
341 else
342 {
343 /* single PDES */
344 if (p.numDim == 3)
345 {
346 Dudley_Assemble_PDE_Single2_3D(p, elements, S, F, A, B, C, D, X, Y);
347 }
348 else if (p.numDim == 2)
349 {
350 Dudley_Assemble_PDE_Single2_2D(p, elements, S, F, A, B, C, D, X, Y);
351 }
352 else
353 {
354 Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
355 }
356 }
357 }
358 else
359 {
360 Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE requires number of equations == number of solutions .");
361 }
362 }
363 blocktimer_increment("Dudley_Assemble_PDE()", blocktimer_start);
364 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26