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

Annotation of /trunk/finley/src/Assemble_PDE.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2548 - (hide annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 17076 byte(s)
Updating copyright notices
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
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 phornby 2078 #include "esysUtils/blocktimer.h"
46 jgs 82 #ifdef _OPENMP
47     #include <omp.h>
48     #endif
49    
50    
51     /**************************************************************/
52    
53 jgs 150 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
54 jgs 82 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
55    
56 gross 798 bool_t reducedIntegrationOrder=FALSE;
57 jgs 150 char error_msg[LenErrorMsg_MAX];
58 gross 798 Assemble_Parameters p;
59 jgs 82 double time0;
60 gross 798 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
61 gross 1028 type_t funcspace;
62 matt 1353 double blocktimer_start = blocktimer_time();
63 gross 798
64 jgs 150 Finley_resetError();
65 jgs 82
66 ksteube 1312 {
67 phornby 1628 #ifdef Paso_MPI
68 ksteube 1312 int iam, numCPUs;
69     iam = elements->elementDistribution->MPIInfo->rank;
70     numCPUs = elements->elementDistribution->MPIInfo->size;
71     #endif
72     }
73    
74 jgs 82 if (nodes==NULL || elements==NULL) return;
75     if (S==NULL && isEmpty(F)) return;
76    
77 matt 1352 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
78 gross 798 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
79     }
80 jgs 82
81 matt 1352 if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
82 gross 798 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
83 jgs 82 }
84 gross 798
85     /* get the functionspace for this assemblage call */
86 gross 1028 funcspace=UNKNOWN;
87 jgs 82 updateFunctionSpaceType(funcspace,A);
88     updateFunctionSpaceType(funcspace,B);
89     updateFunctionSpaceType(funcspace,C);
90     updateFunctionSpaceType(funcspace,D);
91     updateFunctionSpaceType(funcspace,X);
92     updateFunctionSpaceType(funcspace,Y);
93     if (funcspace==UNKNOWN) return; /* all data are empty */
94    
95     /* check if all function spaces are the same */
96     if (! functionSpaceTypeEqual(funcspace,A) ) {
97 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
98 jgs 82 }
99     if (! functionSpaceTypeEqual(funcspace,B) ) {
100 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
101 jgs 82 }
102     if (! functionSpaceTypeEqual(funcspace,C) ) {
103 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
104 jgs 82 }
105     if (! functionSpaceTypeEqual(funcspace,D) ) {
106 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
107 jgs 82 }
108     if (! functionSpaceTypeEqual(funcspace,X) ) {
109 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
110 jgs 82 }
111     if (! functionSpaceTypeEqual(funcspace,Y) ) {
112 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
113 jgs 82 }
114 gross 798 if (! Finley_noError()) return;
115 jgs 82
116     /* check if all function spaces are the same */
117 gross 798 if (funcspace==FINLEY_ELEMENTS) {
118     reducedIntegrationOrder=FALSE;
119     } else if (funcspace==FINLEY_FACE_ELEMENTS) {
120     reducedIntegrationOrder=FALSE;
121     } else if (funcspace==FINLEY_CONTACT_ELEMENTS_1) {
122     reducedIntegrationOrder=FALSE;
123     } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2) {
124     reducedIntegrationOrder=FALSE;
125 gross 1072 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
126     reducedIntegrationOrder=TRUE;
127     } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
128     reducedIntegrationOrder=TRUE;
129     } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
130     reducedIntegrationOrder=TRUE;
131     } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
132     reducedIntegrationOrder=TRUE;
133 gross 798 } else {
134     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");
135     }
136     if (! Finley_noError()) return;
137 jgs 82
138 gross 798 /* set all parameters in p*/
139     Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
140     if (! Finley_noError()) return;
141    
142     /* check if all function spaces are the same */
143    
144 jgs 82 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
145 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
146 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
147 jgs 82 }
148    
149     if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
150 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
151 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
152 jgs 82 }
153    
154     if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
155 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
156 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
157 jgs 82 }
158    
159     if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
160 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
161 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
162 jgs 82 }
163    
164     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
165 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
166 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
167 jgs 82 }
168    
169     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
170 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
171 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
172 jgs 82 }
173    
174     /* check the dimensions: */
175 matt 1352
176 jgs 82 if (p.numEqu==1 && p.numComp==1) {
177     if (!isEmpty(A)) {
178     dimensions[0]=p.numDim;
179     dimensions[1]=p.numDim;
180     if (!isDataPointShapeEqual(A,2,dimensions)) {
181 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
182 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
183 jgs 82 }
184     }
185     if (!isEmpty(B)) {
186     dimensions[0]=p.numDim;
187     if (!isDataPointShapeEqual(B,1,dimensions)) {
188 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
189 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
190 jgs 82 }
191     }
192     if (!isEmpty(C)) {
193     dimensions[0]=p.numDim;
194     if (!isDataPointShapeEqual(C,1,dimensions)) {
195 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
196 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
197 jgs 82 }
198     }
199     if (!isEmpty(D)) {
200     if (!isDataPointShapeEqual(D,0,dimensions)) {
201 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
202 jgs 82 }
203     }
204     if (!isEmpty(X)) {
205     dimensions[0]=p.numDim;
206     if (!isDataPointShapeEqual(X,1,dimensions)) {
207 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
208 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
209 jgs 82 }
210     }
211     if (!isEmpty(Y)) {
212     if (!isDataPointShapeEqual(Y,0,dimensions)) {
213 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
214 jgs 82 }
215 matt 1352 }
216 jgs 82 } else {
217     if (!isEmpty(A)) {
218     dimensions[0]=p.numEqu;
219     dimensions[1]=p.numDim;
220     dimensions[2]=p.numComp;
221     dimensions[3]=p.numDim;
222     if (!isDataPointShapeEqual(A,4,dimensions)) {
223 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
224 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
225 jgs 82 }
226     }
227     if (!isEmpty(B)) {
228     dimensions[0]=p.numEqu;
229     dimensions[1]=p.numDim;
230     dimensions[2]=p.numComp;
231     if (!isDataPointShapeEqual(B,3,dimensions)) {
232 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
233 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
234 jgs 82 }
235     }
236     if (!isEmpty(C)) {
237     dimensions[0]=p.numEqu;
238     dimensions[1]=p.numComp;
239     dimensions[2]=p.numDim;
240     if (!isDataPointShapeEqual(C,3,dimensions)) {
241 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
242 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
243 jgs 82 }
244     }
245     if (!isEmpty(D)) {
246     dimensions[0]=p.numEqu;
247     dimensions[1]=p.numComp;
248     if (!isDataPointShapeEqual(D,2,dimensions)) {
249 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
250 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
251 jgs 82 }
252     }
253     if (!isEmpty(X)) {
254     dimensions[0]=p.numEqu;
255     dimensions[1]=p.numDim;
256     if (!isDataPointShapeEqual(X,2,dimensions)) {
257 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
258 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
259 jgs 82 }
260     }
261     if (!isEmpty(Y)) {
262     dimensions[0]=p.numEqu;
263     if (!isDataPointShapeEqual(Y,1,dimensions)) {
264 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
265 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
266 jgs 82 }
267     }
268     }
269 jgs 150 if (Finley_noError()) {
270 jgs 82 time0=Finley_timer();
271 gross 798 if (p.numEqu == p. numComp) {
272     if (p.numEqu > 1) {
273     /* system of PDESs */
274     if (p.numDim==3) {
275     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
276     Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
277     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
278     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
279     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
280     } else {
281     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
282     }
283     } else {
284     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
285 jgs 82 }
286 gross 798 } else if (p.numDim==2) {
287     if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
288     Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
289     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
290     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
291     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
292     } else {
293     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
294     }
295     } else {
296     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
297     }
298     } else if (p.numDim==2) {
299     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
300     Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);
301     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
302     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
303     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
304     } else {
305     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
306     }
307     } else {
308     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
309     }
310     } else {
311     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
312     }
313     } else {
314     /* single PDES */
315     if (p.numDim==3) {
316     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
317     Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
318     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
319     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
320     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
321     } else {
322     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
323     }
324     } else {
325     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
326     }
327     } else if (p.numDim==2) {
328     if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
329     Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
330     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
331     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
332     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
333     } else {
334     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
335     }
336     } else {
337     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
338     }
339 gross 2156 } else if (p.numDim==1) {
340 gross 798 if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
341     Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);
342     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
343     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
344     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
345     } else {
346     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
347     }
348     } else {
349     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
350     }
351     } else {
352     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
353     }
354     }
355     } else {
356     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions .");
357 jgs 82 }
358 jgs 149 #ifdef Finley_TRACE
359 jgs 82 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
360 jgs 149 #endif
361 jgs 82 }
362 matt 1353 blocktimer_increment("Finley_Assemble_PDE()", blocktimer_start);
363 jgs 82 }
364     /*
365     * $Log$
366 jgs 150 * Revision 1.8 2005/09/15 03:44:21 jgs
367     * Merge of development branch dev-02 back to main trunk on 2005-09-15
368     *
369 jgs 149 * Revision 1.7 2005/09/01 03:31:35 jgs
370     * Merge of development branch dev-02 back to main trunk on 2005-09-01
371     *
372 jgs 147 * Revision 1.6 2005/08/12 01:45:42 jgs
373     * erge of development branch dev-02 back to main trunk on 2005-08-12
374     *
375 jgs 150 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
376     * the solver from finley are put into the standalone package paso now
377     *
378 jgs 149 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
379     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
380     *
381 jgs 147 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
382     * contact element assemblage was called with wrong element table pointer
383     *
384 jgs 123 * Revision 1.5 2005/07/08 04:07:46 jgs
385     * Merge of development branch back to main trunk on 2005-07-08
386     *
387 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
388 jgs 97 * *** empty log message ***
389 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
390     * some changes towards 64 integers in finley
391 jgs 82 *
392 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
393     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
394 jgs 97 *
395 jgs 82 *
396 jgs 123 *
397 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26