/[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 2078 - (show annotations)
Thu Nov 20 16:10:10 2008 UTC (10 years, 11 months ago) by phornby
File MIME type: text/plain
File size: 17076 byte(s)
Two changes.

1. Move blocktimer from escript to esysUtils.
2. Make it possible to link to paso as a DLL or .so.

Should have no effect on 'nix's

In respect of 1., blocktimer had begun to spring up everywhere, so
for the moment I thought it best to move it to the only other library that
pops up all over the place.

In respect of 2., paso needed to be a DLL in order to use the windows intelc /fast
option, which does aggressive multi-file optimisations. Even in its current form, it either
vectorises or parallelises  hundreds more loops in the esys system than appear in the pragmas.

In achieving 2. I have not been too delicate in adding

PASO_DLL_API

declarations to the .h files in paso/src. Only toward the end of the process of
the conversion, when the number of linker errors dropped below 20, say, did I choosy about what
functions in a header I declared PASO_DLL_API. As a result, there are likely to be many routines
declared as external functions symbols that are in fact internal to the paso DLL. 
Why is this an issue?  It prevents the intelc compiler from getting aggressive on the paso module.
With pain there is sometimes gain. At least all the DLL rules in windows give good
(non-microsoft) compiler writers a chance to really shine.

So, if you should see a PASO_DLL_API on a function in a paso header file,
and think to yourself, "that function is only called in paso, why export it?", then feel free to
delete the PASO_DLL_API export declaration.

Here's hoping for no breakage.....
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 "esysUtils/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