/[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 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 9 months ago) by ksteube
File MIME type: text/plain
File size: 17074 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26