/[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 1312 - (hide annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years ago) by ksteube
File MIME type: text/plain
File size: 17114 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26