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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26