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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26