/[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 969 - (show annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 2 months ago) by ksteube
File MIME type: text/plain
File size: 17139 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26