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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4348 - (show annotations)
Tue Apr 2 07:15:43 2013 UTC (6 years ago) by caltinay
File size: 11788 byte(s)
Reapply some of the spelling fixes that got lost....

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/Assemble_PDE.cpp:2682-2741 /branches/pasowrap/dudley/src/Assemble_PDE.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Assemble_PDE.cpp:3871-3891 /branches/restext/finley/src/Assemble_PDE.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Assemble_PDE.cpp:3669-3791 /branches/stage3.0/finley/src/Assemble_PDE.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Assemble_PDE.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Assemble_PDE.cpp:3517-3974 /release/3.0/finley/src/Assemble_PDE.cpp:2591-2601 /trunk/dudley/src/Assemble_PDE.cpp:4257-4344 /trunk/ripley/test/python/dudley/src/Assemble_PDE.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26