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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (9 years, 6 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
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 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26