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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 971 - (hide annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 5 months ago) by ksteube
Original Path: trunk/finley/src/Assemble_PDE.c
File MIME type: text/plain
File size: 16661 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26