/[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 1028 - (hide annotations)
Wed Mar 14 00:15:24 2007 UTC (12 years, 1 month ago) by gross
Original Path: trunk/finley/src/Assemble_PDE.c
File MIME type: text/plain
File size: 16674 byte(s)
modifications to be compliant with _WIN32. The substitutes for asinh, acosh, atanh are still missing (erf will through an exception)
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 gross 798
90 jgs 82 /* check if all function spaces are the same */
91     if (! functionSpaceTypeEqual(funcspace,A) ) {
92 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
93 jgs 82 }
94     if (! functionSpaceTypeEqual(funcspace,B) ) {
95 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
96 jgs 82 }
97     if (! functionSpaceTypeEqual(funcspace,C) ) {
98 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
99 jgs 82 }
100     if (! functionSpaceTypeEqual(funcspace,D) ) {
101 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
102 jgs 82 }
103     if (! functionSpaceTypeEqual(funcspace,X) ) {
104 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
105 jgs 82 }
106     if (! functionSpaceTypeEqual(funcspace,Y) ) {
107 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
108 jgs 82 }
109 gross 798 if (! Finley_noError()) return;
110 jgs 82
111     /* check if all function spaces are the same */
112 gross 798 if (funcspace==FINLEY_ELEMENTS) {
113     reducedIntegrationOrder=FALSE;
114     } else if (funcspace==FINLEY_FACE_ELEMENTS) {
115     reducedIntegrationOrder=FALSE;
116     } else if (funcspace==FINLEY_CONTACT_ELEMENTS_1) {
117     reducedIntegrationOrder=FALSE;
118     } else if (funcspace==FINLEY_CONTACT_ELEMENTS_2) {
119     reducedIntegrationOrder=FALSE;
120     } else {
121     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: assemblage failed because of illegal function space.");
122     }
123     if (! Finley_noError()) return;
124 jgs 82
125 gross 798 /* set all parameters in p*/
126     Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
127     if (! Finley_noError()) return;
128    
129     /* check if all function spaces are the same */
130    
131 jgs 82 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
132 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
133 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
134 jgs 82 }
135    
136     if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
137 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
138 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
139 jgs 82 }
140    
141     if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
142 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
143 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
144 jgs 82 }
145    
146     if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
147 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
148 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
149 jgs 82 }
150    
151     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
152 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
153 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
154 jgs 82 }
155    
156     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
157 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
158 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
159 jgs 82 }
160    
161     /* check the dimensions: */
162    
163     if (p.numEqu==1 && p.numComp==1) {
164     if (!isEmpty(A)) {
165     dimensions[0]=p.numDim;
166     dimensions[1]=p.numDim;
167     if (!isDataPointShapeEqual(A,2,dimensions)) {
168 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
169 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
170 jgs 82 }
171     }
172     if (!isEmpty(B)) {
173     dimensions[0]=p.numDim;
174     if (!isDataPointShapeEqual(B,1,dimensions)) {
175 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
176 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
177 jgs 82 }
178     }
179     if (!isEmpty(C)) {
180     dimensions[0]=p.numDim;
181     if (!isDataPointShapeEqual(C,1,dimensions)) {
182 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
183 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
184 jgs 82 }
185     }
186     if (!isEmpty(D)) {
187     if (!isDataPointShapeEqual(D,0,dimensions)) {
188 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
189 jgs 82 }
190     }
191     if (!isEmpty(X)) {
192     dimensions[0]=p.numDim;
193     if (!isDataPointShapeEqual(X,1,dimensions)) {
194 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
195 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
196 jgs 82 }
197     }
198     if (!isEmpty(Y)) {
199     if (!isDataPointShapeEqual(Y,0,dimensions)) {
200 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
201 jgs 82 }
202     }
203     } else {
204     if (!isEmpty(A)) {
205     dimensions[0]=p.numEqu;
206     dimensions[1]=p.numDim;
207     dimensions[2]=p.numComp;
208     dimensions[3]=p.numDim;
209     if (!isDataPointShapeEqual(A,4,dimensions)) {
210 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
211 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
212 jgs 82 }
213     }
214     if (!isEmpty(B)) {
215     dimensions[0]=p.numEqu;
216     dimensions[1]=p.numDim;
217     dimensions[2]=p.numComp;
218     if (!isDataPointShapeEqual(B,3,dimensions)) {
219 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
220 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
221 jgs 82 }
222     }
223     if (!isEmpty(C)) {
224     dimensions[0]=p.numEqu;
225     dimensions[1]=p.numComp;
226     dimensions[2]=p.numDim;
227     if (!isDataPointShapeEqual(C,3,dimensions)) {
228 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
229 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
230 jgs 82 }
231     }
232     if (!isEmpty(D)) {
233     dimensions[0]=p.numEqu;
234     dimensions[1]=p.numComp;
235     if (!isDataPointShapeEqual(D,2,dimensions)) {
236 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
237 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
238 jgs 82 }
239     }
240     if (!isEmpty(X)) {
241     dimensions[0]=p.numEqu;
242     dimensions[1]=p.numDim;
243     if (!isDataPointShapeEqual(X,2,dimensions)) {
244 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
245 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
246 jgs 82 }
247     }
248     if (!isEmpty(Y)) {
249     dimensions[0]=p.numEqu;
250     if (!isDataPointShapeEqual(Y,1,dimensions)) {
251 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
252 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
253 jgs 82 }
254     }
255     }
256 jgs 150 if (Finley_noError()) {
257 jgs 82 time0=Finley_timer();
258 gross 798 if (p.numEqu == p. numComp) {
259     if (p.numEqu > 1) {
260     /* system of PDESs */
261     if (p.numDim==3) {
262     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
263     Finley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
264     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
265     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
266     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
267     } else {
268     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
269     }
270     } else {
271     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
272 jgs 82 }
273 gross 798 } else if (p.numDim==2) {
274     if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
275     Finley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
276     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
277     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
278     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
279     } else {
280     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
281     }
282     } else {
283     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
284     }
285     } else if (p.numDim==2) {
286     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
287     Finley_Assemble_PDE_System2_1D(p,elements,S,F,A,B,C,D,X,Y);
288     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
289     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
290     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
291     } else {
292     Finley_Assemble_PDE_System2_C(p,elements,S,F,D,Y);
293     }
294     } else {
295     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
296     }
297     } else {
298     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
299     }
300     } else {
301     /* single PDES */
302     if (p.numDim==3) {
303     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
304     Finley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
305     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
306     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
307     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
308     } else {
309     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
310     }
311     } else {
312     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
313     }
314     } else if (p.numDim==2) {
315     if ((p.row_NS == p.col_NS) && (p.row_NS == p.row_NN) && (p.col_NS == p.col_NN )) {
316     Finley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
317     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
318     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
319     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
320     } else {
321     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
322     }
323     } else {
324     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
325     }
326     } else if (p.numDim==2) {
327     if (p.row_NS == p.col_NS && p.row_NS == p.row_NN && p.col_NS == p.col_NN ) {
328     Finley_Assemble_PDE_Single2_1D(p,elements,S,F,A,B,C,D,X,Y);
329     } else if ( p.row_NS == p.col_NS && 2*p.row_NS == p.row_NN && 2*p.col_NS == p.col_NN ) {
330     if ( !isEmpty(A) || !isEmpty(B) || !isEmpty(C) || !isEmpty(X) ) {
331     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: Contact elements require A, B, C and X to be empty.");
332     } else {
333     Finley_Assemble_PDE_Single2_C(p,elements,S,F,D,Y);
334     }
335     } else {
336     Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE supports numShape=NumNodes or 2*numShape=NumNodes only.");
337     }
338     } else {
339     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE supports spatial dimensions 1,2,3 only.");
340     }
341     }
342     } else {
343     Finley_setError(VALUE_ERROR,"Finley_Assemble_PDE requires number of equations == number of solutions .");
344 jgs 82 }
345 jgs 149 #ifdef Finley_TRACE
346 jgs 82 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
347 jgs 149 #endif
348 jgs 82 }
349     }
350     /*
351     * $Log$
352 jgs 150 * Revision 1.8 2005/09/15 03:44:21 jgs
353     * Merge of development branch dev-02 back to main trunk on 2005-09-15
354     *
355 jgs 149 * Revision 1.7 2005/09/01 03:31:35 jgs
356     * Merge of development branch dev-02 back to main trunk on 2005-09-01
357     *
358 jgs 147 * Revision 1.6 2005/08/12 01:45:42 jgs
359     * erge of development branch dev-02 back to main trunk on 2005-08-12
360     *
361 jgs 150 * Revision 1.5.2.3 2005/09/07 06:26:17 gross
362     * the solver from finley are put into the standalone package paso now
363     *
364 jgs 149 * Revision 1.5.2.2 2005/08/24 02:02:18 gross
365     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
366     *
367 jgs 147 * Revision 1.5.2.1 2005/08/03 08:54:27 gross
368     * contact element assemblage was called with wrong element table pointer
369     *
370 jgs 123 * Revision 1.5 2005/07/08 04:07:46 jgs
371     * Merge of development branch back to main trunk on 2005-07-08
372     *
373 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
374 jgs 97 * *** empty log message ***
375 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
376     * some changes towards 64 integers in finley
377 jgs 82 *
378 jgs 123 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
379     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
380 jgs 97 *
381 jgs 82 *
382 jgs 123 *
383 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26