/[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 1072 - (hide annotations)
Thu Mar 29 06:44:30 2007 UTC (12 years, 4 months ago) by gross
Original Path: trunk/finley/src/Assemble_PDE.c
File MIME type: text/plain
File size: 17053 byte(s)
PDE assemblage for reduced integration order + tests added.


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 gross 1028 type_t funcspace;
65 gross 798
66 jgs 150 Finley_resetError();
67 jgs 82
68     if (nodes==NULL || elements==NULL) return;
69     if (S==NULL && isEmpty(F)) return;
70    
71 gross 798 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
72     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
73     }
74 jgs 82
75 gross 798 if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
76     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
77 jgs 82 }
78 gross 798
79     /* get the functionspace for this assemblage call */
80 gross 1028 funcspace=UNKNOWN;
81 jgs 82 updateFunctionSpaceType(funcspace,A);
82     updateFunctionSpaceType(funcspace,B);
83     updateFunctionSpaceType(funcspace,C);
84     updateFunctionSpaceType(funcspace,D);
85     updateFunctionSpaceType(funcspace,X);
86     updateFunctionSpaceType(funcspace,Y);
87     if (funcspace==UNKNOWN) return; /* all data are empty */
88    
89     /* 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 gross 1072 } else if (funcspace==FINLEY_REDUCED_ELEMENTS) {
120     reducedIntegrationOrder=TRUE;
121     } else if (funcspace==FINLEY_REDUCED_FACE_ELEMENTS) {
122     reducedIntegrationOrder=TRUE;
123     } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_1) {
124     reducedIntegrationOrder=TRUE;
125     } else if (funcspace==FINLEY_REDUCED_CONTACT_ELEMENTS_2) {
126     reducedIntegrationOrder=TRUE;
127 gross 798 } else {
128     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");
129     }
130     if (! Finley_noError()) return;
131 jgs 82
132 gross 798 /* set all parameters in p*/
133     Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
134     if (! Finley_noError()) return;
135    
136     /* check if all function spaces are the same */
137    
138 jgs 82 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
139 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
140 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
141 jgs 82 }
142    
143     if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
144 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
145 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
146 jgs 82 }
147    
148     if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
149 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
150 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
151 jgs 82 }
152    
153     if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
154 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
155 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
156 jgs 82 }
157    
158     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
159 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
160 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
161 jgs 82 }
162    
163     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
164 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
165 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
166 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
176 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
177 jgs 82 }
178     }
179     if (!isEmpty(B)) {
180     dimensions[0]=p.numDim;
181     if (!isDataPointShapeEqual(B,1,dimensions)) {
182 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
183 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
184 jgs 82 }
185     }
186     if (!isEmpty(C)) {
187     dimensions[0]=p.numDim;
188     if (!isDataPointShapeEqual(C,1,dimensions)) {
189 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
190 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
191 jgs 82 }
192     }
193     if (!isEmpty(D)) {
194     if (!isDataPointShapeEqual(D,0,dimensions)) {
195 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
196 jgs 82 }
197     }
198     if (!isEmpty(X)) {
199     dimensions[0]=p.numDim;
200     if (!isDataPointShapeEqual(X,1,dimensions)) {
201 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
202 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
203 jgs 82 }
204     }
205     if (!isEmpty(Y)) {
206     if (!isDataPointShapeEqual(Y,0,dimensions)) {
207 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
208 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
218 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
219 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
227 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
228 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
236 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
237 jgs 82 }
238     }
239     if (!isEmpty(D)) {
240     dimensions[0]=p.numEqu;
241     dimensions[1]=p.numComp;
242     if (!isDataPointShapeEqual(D,2,dimensions)) {
243 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
244 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
245 jgs 82 }
246     }
247     if (!isEmpty(X)) {
248     dimensions[0]=p.numEqu;
249     dimensions[1]=p.numDim;
250     if (!isDataPointShapeEqual(X,2,dimensions)) {
251 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
252 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
253 jgs 82 }
254     }
255     if (!isEmpty(Y)) {
256     dimensions[0]=p.numEqu;
257     if (!isDataPointShapeEqual(Y,1,dimensions)) {
258 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
259 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
260 jgs 82 }
261     }
262     }
263 jgs 150 if (Finley_noError()) {
264 jgs 82 time0=Finley_timer();
265 gross 798 if (p.numEqu == p. numComp) {
266     if (p.numEqu > 1) {
267     /* system of PDESs */
268     if (p.numDim==3) {
269     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
270     Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
271     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
272     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
273     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
274     } else {
275     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
276     }
277     } else {
278     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
279 jgs 82 }
280 gross 798 } else if (p.numDim==2) {
281     if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
282     Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
283     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
284     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
285     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
286     } else {
287     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
288     }
289     } else {
290     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
291     }
292     } else if (p.numDim==2) {
293     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
294     Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);
295     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
296     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
297     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
298     } else {
299     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
300     }
301     } else {
302     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
303     }
304     } else {
305     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
306     }
307     } else {
308     /* single PDES */
309     if (p.numDim==3) {
310     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
311     Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
312     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
313     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
314     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
315     } else {
316     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
317     }
318     } else {
319     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
320     }
321     } else if (p.numDim==2) {
322     if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
323     Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
324     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
325     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
326     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
327     } else {
328     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
329     }
330     } else {
331     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
332     }
333     } else if (p.numDim==2) {
334     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
335     Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);
336     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
337     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
338     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
339     } else {
340     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
341     }
342     } else {
343     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
344     }
345     } else {
346     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
347     }
348     }
349     } else {
350     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions .");
351 jgs 82 }
352 jgs 149 #ifdef Finley_TRACE
353 jgs 82 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
354 jgs 149 #endif
355 jgs 82 }
356     }
357     /*
358     * $Log$
359 jgs 150 * Revision 1.8 2005/09/15 03:44:21 jgs
360     * Merge of development branch dev-02 back to main trunk on 2005-09-15
361     *
362 jgs 149 * Revision 1.7 2005/09/01 03:31:35 jgs
363     * Merge of development branch dev-02 back to main trunk on 2005-09-01
364     *
365 jgs 147 * Revision 1.6 2005/08/12 01:45:42 jgs
366     * erge of development branch dev-02 back to main trunk on 2005-08-12
367     *
368 jgs 150 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
369     * the solver from finley are put into the standalone package paso now
370     *
371 jgs 149 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
372     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
373     *
374 jgs 147 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
375     * contact element assemblage was called with wrong element table pointer
376     *
377 jgs 123 * Revision 1.5 2005/07/08 04:07:46 jgs
378     * Merge of development branch back to main trunk on 2005-07-08
379     *
380 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
381 jgs 97 * *** empty log message ***
382 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
383     * some changes towards 64 integers in finley
384 jgs 82 *
385 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
386     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
387 jgs 97 *
388 jgs 82 *
389 jgs 123 *
390 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26