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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 971 - (show annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 8 months ago) by ksteube
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 /*
2 ************************************************************
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 */
12
13
14 /**************************************************************/
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 /* Author: gross@access.edu.au */
43 /* Version: $Id$ */
44
45 /**************************************************************/
46
47 #include "Assemble.h"
48 #include "Util.h"
49 #ifdef _OPENMP
50 #include <omp.h>
51 #endif
52
53
54 /**************************************************************/
55
56 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
57 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
58
59 bool_t reducedIntegrationOrder=FALSE;
60 char error_msg[LenErrorMsg_MAX];
61 Assemble_Parameters p;
62 double time0;
63 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
64
65 Finley_resetError();
66
67 if (nodes==NULL || elements==NULL) return;
68 if (S==NULL && isEmpty(F)) return;
69
70 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
74 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 }
77
78 /* get the functionspace for this assemblage call */
79 type_t funcspace=UNKNOWN;
80 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
89 /* check if all function spaces are the same */
90 if (! functionSpaceTypeEqual(funcspace,A) ) {
91 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
92 }
93 if (! functionSpaceTypeEqual(funcspace,B) ) {
94 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
95 }
96 if (! functionSpaceTypeEqual(funcspace,C) ) {
97 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
98 }
99 if (! functionSpaceTypeEqual(funcspace,D) ) {
100 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
101 }
102 if (! functionSpaceTypeEqual(funcspace,X) ) {
103 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
104 }
105 if (! functionSpaceTypeEqual(funcspace,Y) ) {
106 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
107 }
108 if (! Finley_noError()) return;
109
110 /* check if all function spaces are the same */
111 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
124 /* 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 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
131 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
132 Finley_setError(TYPE_ERROR,error_msg);
133 }
134
135 if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
136 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
137 Finley_setError(TYPE_ERROR,error_msg);
138 }
139
140 if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
141 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
142 Finley_setError(TYPE_ERROR,error_msg);
143 }
144
145 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
146 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
147 Finley_setError(TYPE_ERROR,error_msg);
148 }
149
150 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
151 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
152 Finley_setError(TYPE_ERROR,error_msg);
153 }
154
155 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
156 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
157 Finley_setError(TYPE_ERROR,error_msg);
158 }
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 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
168 Finley_setError(TYPE_ERROR,error_msg);
169 }
170 }
171 if (!isEmpty(B)) {
172 dimensions[0]=p.numDim;
173 if (!isDataPointShapeEqual(B,1,dimensions)) {
174 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
175 Finley_setError(TYPE_ERROR,error_msg);
176 }
177 }
178 if (!isEmpty(C)) {
179 dimensions[0]=p.numDim;
180 if (!isDataPointShapeEqual(C,1,dimensions)) {
181 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
182 Finley_setError(TYPE_ERROR,error_msg);
183 }
184 }
185 if (!isEmpty(D)) {
186 if (!isDataPointShapeEqual(D,0,dimensions)) {
187 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
188 }
189 }
190 if (!isEmpty(X)) {
191 dimensions[0]=p.numDim;
192 if (!isDataPointShapeEqual(X,1,dimensions)) {
193 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
194 Finley_setError(TYPE_ERROR,error_msg);
195 }
196 }
197 if (!isEmpty(Y)) {
198 if (!isDataPointShapeEqual(Y,0,dimensions)) {
199 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
200 }
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 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
210 Finley_setError(TYPE_ERROR,error_msg);
211 }
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 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
219 Finley_setError(TYPE_ERROR,error_msg);
220 }
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 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
228 Finley_setError(TYPE_ERROR,error_msg);
229 }
230 }
231 if (!isEmpty(D)) {
232 dimensions[0]=p.numEqu;
233 dimensions[1]=p.numComp;
234 if (!isDataPointShapeEqual(D,2,dimensions)) {
235 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
236 Finley_setError(TYPE_ERROR,error_msg);
237 }
238 }
239 if (!isEmpty(X)) {
240 dimensions[0]=p.numEqu;
241 dimensions[1]=p.numDim;
242 if (!isDataPointShapeEqual(X,2,dimensions)) {
243 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
244 Finley_setError(TYPE_ERROR,error_msg);
245 }
246 }
247 if (!isEmpty(Y)) {
248 dimensions[0]=p.numEqu;
249 if (!isDataPointShapeEqual(Y,1,dimensions)) {
250 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
251 Finley_setError(TYPE_ERROR,error_msg);
252 }
253 }
254 }
255 if (Finley_noError()) {
256 time0=Finley_timer();
257 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 }
272 } 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 }
344 #ifdef Finley_TRACE
345 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
346 #endif
347 }
348 }
349 /*
350 * $Log$
351 * 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 * 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 * 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 * 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 * 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 * 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 * 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 * Revision 1.4 2004/12/15 07:08:32 jgs
373 * *** empty log message ***
374 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
375 * some changes towards 64 integers in finley
376 *
377 * 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 *
380 *
381 *
382 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26