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

Contents of /trunk/finley/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, 3 months ago) by gross
File MIME type: text/plain
File size: 17251 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
17 /* assembles the system of numEq PDEs into the stiffness matrix S and right hand side F */
18
19 /* -div(A*grad u)-div(B*u)+C*grad u + D*u= -div X + Y */
20
21 /* -(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 */
22
23 /* u has numComp components. */
24
25 /* Shape of the coefficients: */
26
27 /* A = numEqu x numDim x numComp x numDim */
28 /* B = numDim x numEqu x numComp */
29 /* C = numEqu x numDim x numComp */
30 /* D = numEqu x numComp */
31 /* X = numEqu x numDim */
32 /* Y = numEqu */
33
34 /* The coefficients A,B,C,D,X and Y have to be defined on the integartion points or not present (=NULL). */
35
36 /* S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
37 /* the right hand side of the PDE is not processed. */
38
39 /* The routine does not consider any boundary conditions. */
40
41 /**************************************************************/
42
43 #include "Assemble.h"
44 #include "Util.h"
45 #include "esysUtils/blocktimer.h"
46 #ifdef _OPENMP
47 #include <omp.h>
48 #endif
49
50
51 /**************************************************************/
52
53 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
54 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
55
56 bool_t reducedIntegrationOrder=FALSE;
57 char error_msg[LenErrorMsg_MAX];
58 Finley_Assemble_Parameters p;
59 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
60 type_t funcspace;
61 double blocktimer_start = blocktimer_time();
62
63 Finley_resetError();
64
65 if (nodes==NULL || elements==NULL) return;
66 if (S==NULL && isEmpty(F)) return;
67
68 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
69 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
70 }
71
72 if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
73 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
74 }
75
76 /* get the functionspace for this assemblage call */
77 funcspace=UNKNOWN;
78 updateFunctionSpaceType(funcspace,A);
79 updateFunctionSpaceType(funcspace,B);
80 updateFunctionSpaceType(funcspace,C);
81 updateFunctionSpaceType(funcspace,D);
82 updateFunctionSpaceType(funcspace,X);
83 updateFunctionSpaceType(funcspace,Y);
84 if (funcspace==UNKNOWN) return; /* all data are empty */
85
86 /* check if all function spaces are the same */
87 if (! functionSpaceTypeEqual(funcspace,A) ) {
88 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
89 }
90 if (! functionSpaceTypeEqual(funcspace,B) ) {
91 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
92 }
93 if (! functionSpaceTypeEqual(funcspace,C) ) {
94 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
95 }
96 if (! functionSpaceTypeEqual(funcspace,D) ) {
97 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
98 }
99 if (! functionSpaceTypeEqual(funcspace,X) ) {
100 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
101 }
102 if (! functionSpaceTypeEqual(funcspace,Y) ) {
103 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
104 }
105 if (! Finley_noError()) return;
106
107 /* check if all function spaces are the same */
108 if (funcspace==FINLEY_ELEMENTS) {
109 reducedIntegrationOrder=FALSE;
110 } else if (funcspace==FINLEY_FACE_ELEMENTS) {
111 reducedIntegrationOrder=FALSE;
112 } else if (funcspace==FINLEY_CONTACT_ELEMENTS_1) {
113 reducedIntegrationOrder=FALSE;
114 } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2) {
115 reducedIntegrationOrder=FALSE;
116 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
117 reducedIntegrationOrder=TRUE;
118 } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
119 reducedIntegrationOrder=TRUE;
120 } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
121 reducedIntegrationOrder=TRUE;
122 } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
123 reducedIntegrationOrder=TRUE;
124 } else if (funcspace==FINLEY_POINTS) {
125 reducedIntegrationOrder=FALSE;
126 } else {
127 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");
128 }
129 if (! Finley_noError()) return;
130
131 /* set all parameters in p*/
132 Finley_Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
133 if (! Finley_noError()) return;
134
135 /* check if all function spaces are the same */
136
137 if (! numSamplesEqual(A,p.numQuadTotal,elements->numElements) ) {
138 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuadTotal,elements->numElements);
139 Finley_setError(TYPE_ERROR,error_msg);
140 }
141
142 if (! numSamplesEqual(B,p.numQuadTotal,elements->numElements) ) {
143 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuadTotal,elements->numElements);
144 Finley_setError(TYPE_ERROR,error_msg);
145 }
146
147 if (! numSamplesEqual(C,p.numQuadTotal,elements->numElements) ) {
148 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuadTotal,elements->numElements);
149 Finley_setError(TYPE_ERROR,error_msg);
150 }
151
152 if (! numSamplesEqual(D,p.numQuadTotal,elements->numElements) ) {
153 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuadTotal,elements->numElements);
154 Finley_setError(TYPE_ERROR,error_msg);
155 }
156
157 if (! numSamplesEqual(X,p.numQuadTotal,elements->numElements) ) {
158 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuadTotal,elements->numElements);
159 Finley_setError(TYPE_ERROR,error_msg);
160 }
161
162 if (! numSamplesEqual(Y,p.numQuadTotal,elements->numElements) ) {
163 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuadTotal,elements->numElements);
164 Finley_setError(TYPE_ERROR,error_msg);
165 }
166
167 /* check the dimensions: */
168
169 if (p.numEqu==1 && p.numComp==1) {
170 if (!isEmpty(A)) {
171 dimensions[0]=p.numDim;
172 dimensions[1]=p.numDim;
173 if (!isDataPointShapeEqual(A,2,dimensions)) {
174 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
175 Finley_setError(TYPE_ERROR,error_msg);
176 }
177 }
178 if (!isEmpty(B)) {
179 dimensions[0]=p.numDim;
180 if (!isDataPointShapeEqual(B,1,dimensions)) {
181 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
182 Finley_setError(TYPE_ERROR,error_msg);
183 }
184 }
185 if (!isEmpty(C)) {
186 dimensions[0]=p.numDim;
187 if (!isDataPointShapeEqual(C,1,dimensions)) {
188 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
189 Finley_setError(TYPE_ERROR,error_msg);
190 }
191 }
192 if (!isEmpty(D)) {
193 if (!isDataPointShapeEqual(D,0,dimensions)) {
194 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
195 }
196 }
197 if (!isEmpty(X)) {
198 dimensions[0]=p.numDim;
199 if (!isDataPointShapeEqual(X,1,dimensions)) {
200 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
201 Finley_setError(TYPE_ERROR,error_msg);
202 }
203 }
204 if (!isEmpty(Y)) {
205 if (!isDataPointShapeEqual(Y,0,dimensions)) {
206 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
207 }
208 }
209 } else {
210 if (!isEmpty(A)) {
211 dimensions[0]=p.numEqu;
212 dimensions[1]=p.numDim;
213 dimensions[2]=p.numComp;
214 dimensions[3]=p.numDim;
215 if (!isDataPointShapeEqual(A,4,dimensions)) {
216 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
217 Finley_setError(TYPE_ERROR,error_msg);
218 }
219 }
220 if (!isEmpty(B)) {
221 dimensions[0]=p.numEqu;
222 dimensions[1]=p.numDim;
223 dimensions[2]=p.numComp;
224 if (!isDataPointShapeEqual(B,3,dimensions)) {
225 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
226 Finley_setError(TYPE_ERROR,error_msg);
227 }
228 }
229 if (!isEmpty(C)) {
230 dimensions[0]=p.numEqu;
231 dimensions[1]=p.numComp;
232 dimensions[2]=p.numDim;
233 if (!isDataPointShapeEqual(C,3,dimensions)) {
234 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
235 Finley_setError(TYPE_ERROR,error_msg);
236 }
237 }
238 if (!isEmpty(D)) {
239 dimensions[0]=p.numEqu;
240 dimensions[1]=p.numComp;
241 if (!isDataPointShapeEqual(D,2,dimensions)) {
242 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
243 Finley_setError(TYPE_ERROR,error_msg);
244 }
245 }
246 if (!isEmpty(X)) {
247 dimensions[0]=p.numEqu;
248 dimensions[1]=p.numDim;
249 if (!isDataPointShapeEqual(X,2,dimensions)) {
250 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
251 Finley_setError(TYPE_ERROR,error_msg);
252 }
253 }
254 if (!isEmpty(Y)) {
255 dimensions[0]=p.numEqu;
256 if (!isDataPointShapeEqual(Y,1,dimensions)) {
257 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
258 Finley_setError(TYPE_ERROR,error_msg);
259 }
260 }
261 }
262 if (Finley_noError()) {
263 if (p.numEqu == p. numComp) {
264 if (p.numEqu > 1) {
265 /* system of PDESs */
266 if (p.numDim==3) {
267 if ( p.numSides == 1 ) {
268
269 if (funcspace==FINLEY_POINTS) {
270 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
271 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
272 } else {
273 Finley_Assemble_PDE_Points(p, elements,S,F, D, Y);
274 }
275 } else {
276 Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
277 }
278
279 } else if ( p.numSides == 2 ) {
280 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
281 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
282 } else {
283 Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
284 }
285 } else {
286 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
287 }
288 } else if (p.numDim==2) {
289 if ( p.numSides == 1 ) {
290 if (funcspace==FINLEY_POINTS) {
291 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
292 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
293 } else {
294 Finley_Assemble_PDE_Points(p, elements,S,F, D, Y);
295 }
296 } else {
297 Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
298 }
299 } else if ( p.numSides == 2 ) {
300 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
301 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
302 } else {
303 Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
304 }
305 } else {
306 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
307 }
308 } else if (p.numDim==1) {
309 if ( p.numSides == 1 ) {
310 if (funcspace==FINLEY_POINTS) {
311 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
312 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
313 } else {
314 Finley_Assemble_PDE_Points(p, elements,S,F, D, Y);
315 }
316 } else {
317 Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);
318 }
319 } else if ( p.numSides == 2 ) {
320 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
321 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
322 } else {
323 Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
324 }
325 } else {
326 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
327 }
328 } else {
329 Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
330 }
331 } else {
332 /* single PDES */
333 if (p.numDim==3) {
334 if ( p.numSides == 1 ) {
335 if (funcspace==FINLEY_POINTS) {
336 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
337 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
338 } else {
339 Finley_Assemble_PDE_Points(p, elements,S,F, D, Y);
340 }
341 } else {
342 Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
343 }
344 } else if ( p.numSides == 2 ) {
345 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
346 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
347 } else {
348 Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
349 }
350 } else {
351 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
352 }
353 } else if (p.numDim==2) {
354 if ( p.numSides == 1 ) {
355 if (funcspace==FINLEY_POINTS) {
356 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
357 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
358 } else {
359 Finley_Assemble_PDE_Points(p, elements,S,F, D, Y);
360 }
361 } else {
362 Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
363 }
364 } else if ( p.numSides == 2 ) {
365 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
366 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
367 } else {
368 Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
369 }
370 } else {
371 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
372 }
373 } else if (p.numDim==1) {
374 if ( p.numSides == 1 ) {
375 if (funcspace==FINLEY_POINTS) {
376 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
377 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Point elements require A, B, C and X to be empty.");
378 } else {
379 Finley_Assemble_PDE_Points(p, elements,S,F, D, Y);
380 }
381 } else {
382 Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);
383 }
384 } else if ( p.numSides == 2 ) {
385 if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
386 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
387 } else {
388 Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
389 }
390 } else {
391 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
392 }
393 } else {
394 Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
395 }
396 }
397 } else {
398 Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions .");
399 }
400 }
401 blocktimer_increment("Finley_Assemble_PDE()", blocktimer_start);
402 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26