/[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 3522 - (show annotations)
Tue May 24 00:57:58 2011 UTC (8 years, 4 months ago) by gross
File MIME type: text/plain
File size: 11598 byte(s)
(almost) full support for Point elements
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 if (funcspace == DUDLEY_POINTS)
138 {
139 reducedIntegrationOrder = TRUE;
140 }
141 else
142 {
143 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: assemblage failed because of illegal function space.");
144 }
145 if (!Dudley_noError())
146 return;
147
148 /* set all parameters in p */
149 Dudley_Assemble_getAssembleParameters(nodes, elements, S, F, reducedIntegrationOrder, &p);
150 if (!Dudley_noError())
151 return;
152
153 /* check if all function spaces are the same */
154
155 if (!numSamplesEqual(A, p.numQuad, elements->numElements))
156 {
157 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)", p.numQuad,
158 elements->numElements);
159 Dudley_setError(TYPE_ERROR, error_msg);
160 }
161
162 if (!numSamplesEqual(B, p.numQuad, elements->numElements))
163 {
164 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)", p.numQuad,
165 elements->numElements);
166 Dudley_setError(TYPE_ERROR, error_msg);
167 }
168
169 if (!numSamplesEqual(C, p.numQuad, elements->numElements))
170 {
171 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)", p.numQuad,
172 elements->numElements);
173 Dudley_setError(TYPE_ERROR, error_msg);
174 }
175
176 if (!numSamplesEqual(D, p.numQuad, elements->numElements))
177 {
178 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)", p.numQuad,
179 elements->numElements);
180 Dudley_setError(TYPE_ERROR, error_msg);
181 }
182
183 if (!numSamplesEqual(X, p.numQuad, elements->numElements))
184 {
185 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)", p.numQuad,
186 elements->numElements);
187 Dudley_setError(TYPE_ERROR, error_msg);
188 }
189
190 if (!numSamplesEqual(Y, p.numQuad, elements->numElements))
191 {
192 sprintf(error_msg, "Dudley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)", p.numQuad,
193 elements->numElements);
194 Dudley_setError(TYPE_ERROR, error_msg);
195 }
196
197 /* check the dimensions: */
198
199 if (p.numEqu == 1 && p.numComp == 1)
200 {
201 if (!isEmpty(A))
202 {
203 dimensions[0] = p.numDim;
204 dimensions[1] = p.numDim;
205 if (!isDataPointShapeEqual(A, 2, dimensions))
206 {
207 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",
208 dimensions[0], dimensions[1]);
209 Dudley_setError(TYPE_ERROR, error_msg);
210 }
211 }
212 if (!isEmpty(B))
213 {
214 dimensions[0] = p.numDim;
215 if (!isDataPointShapeEqual(B, 1, dimensions))
216 {
217 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient B: illegal shape (%d,)", dimensions[0]);
218 Dudley_setError(TYPE_ERROR, error_msg);
219 }
220 }
221 if (!isEmpty(C))
222 {
223 dimensions[0] = p.numDim;
224 if (!isDataPointShapeEqual(C, 1, dimensions))
225 {
226 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient C, expected shape (%d,)", dimensions[0]);
227 Dudley_setError(TYPE_ERROR, error_msg);
228 }
229 }
230 if (!isEmpty(D))
231 {
232 if (!isDataPointShapeEqual(D, 0, dimensions))
233 {
234 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficient D, rank 0 expected.");
235 }
236 }
237 if (!isEmpty(X))
238 {
239 dimensions[0] = p.numDim;
240 if (!isDataPointShapeEqual(X, 1, dimensions))
241 {
242 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient X, expected shape (%d,", dimensions[0]);
243 Dudley_setError(TYPE_ERROR, error_msg);
244 }
245 }
246 if (!isEmpty(Y))
247 {
248 if (!isDataPointShapeEqual(Y, 0, dimensions))
249 {
250 Dudley_setError(TYPE_ERROR, "Dudley_Assemble_PDE: coefficient Y, rank 0 expected.");
251 }
252 }
253 }
254 else
255 {
256 if (!isEmpty(A))
257 {
258 dimensions[0] = p.numEqu;
259 dimensions[1] = p.numDim;
260 dimensions[2] = p.numComp;
261 dimensions[3] = p.numDim;
262 if (!isDataPointShapeEqual(A, 4, dimensions))
263 {
264 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)", dimensions[0],
265 dimensions[1], dimensions[2], dimensions[3]);
266 Dudley_setError(TYPE_ERROR, error_msg);
267 }
268 }
269 if (!isEmpty(B))
270 {
271 dimensions[0] = p.numEqu;
272 dimensions[1] = p.numDim;
273 dimensions[2] = p.numComp;
274 if (!isDataPointShapeEqual(B, 3, dimensions))
275 {
276 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)", dimensions[0],
277 dimensions[1], dimensions[2]);
278 Dudley_setError(TYPE_ERROR, error_msg);
279 }
280 }
281 if (!isEmpty(C))
282 {
283 dimensions[0] = p.numEqu;
284 dimensions[1] = p.numComp;
285 dimensions[2] = p.numDim;
286 if (!isDataPointShapeEqual(C, 3, dimensions))
287 {
288 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)", dimensions[0],
289 dimensions[1], dimensions[2]);
290 Dudley_setError(TYPE_ERROR, error_msg);
291 }
292 }
293 if (!isEmpty(D))
294 {
295 dimensions[0] = p.numEqu;
296 dimensions[1] = p.numComp;
297 if (!isDataPointShapeEqual(D, 2, dimensions))
298 {
299 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient D, expected shape (%d,%d)", dimensions[0],
300 dimensions[1]);
301 Dudley_setError(TYPE_ERROR, error_msg);
302 }
303 }
304 if (!isEmpty(X))
305 {
306 dimensions[0] = p.numEqu;
307 dimensions[1] = p.numDim;
308 if (!isDataPointShapeEqual(X, 2, dimensions))
309 {
310 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient X, expected shape (%d,%d)", dimensions[0],
311 dimensions[1]);
312 Dudley_setError(TYPE_ERROR, error_msg);
313 }
314 }
315 if (!isEmpty(Y))
316 {
317 dimensions[0] = p.numEqu;
318 if (!isDataPointShapeEqual(Y, 1, dimensions))
319 {
320 sprintf(error_msg, "Dudley_Assemble_PDE: coefficient Y, expected shape (%d,)", dimensions[0]);
321 Dudley_setError(TYPE_ERROR, error_msg);
322 }
323 }
324 }
325 if (Dudley_noError())
326 {
327 if (funcspace==DUDLEY_POINTS) {
328 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
329 Dudley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
330 } else {
331 Dudley_Assemble_PDE_Points(p, elements,S,F, D, Y);
332 }
333 }
334 else
335 {
336 if (p.numEqu == p.numComp)
337 {
338 if (p.numEqu > 1)
339 {
340 /* system of PDESs */
341 if (p.numDim == 3)
342 {
343 Dudley_Assemble_PDE_System2_3D(p, elements, S, F, A, B, C, D, X, Y);
344 }
345 else if (p.numDim == 2)
346 {
347 Dudley_Assemble_PDE_System2_2D(p, elements, S, F, A, B, C, D, X, Y);
348 }
349 else
350 {
351 Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
352 }
353 }
354 else
355 {
356 /* single PDES */
357 if (p.numDim == 3)
358 {
359 Dudley_Assemble_PDE_Single2_3D(p, elements, S, F, A, B, C, D, X, Y);
360 }
361 else if (p.numDim == 2)
362 {
363 Dudley_Assemble_PDE_Single2_2D(p, elements, S, F, A, B, C, D, X, Y);
364 }
365 else
366 {
367 Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
368 }
369 }
370 }
371 else
372 {
373 Dudley_setError(VALUE_ERROR, "Dudley_Assemble_PDE requires number of equations == number of solutions .");
374 }
375 }
376 }
377 blocktimer_increment("Dudley_Assemble_PDE()", blocktimer_start);
378 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26