/[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 1385 - (show annotations)
Fri Jan 11 07:33:30 2008 UTC (13 years, 1 month ago) by trankine
File MIME type: text/plain
File size: 16963 byte(s)
This is a commit to the trunk of the current windows version undergoing debugging.
This trunk will shortly be moved to the branches area
of the repository, and the saved trunk version currently over in branches restored to the trunk.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26