/[escript]/trunk/finley/src/Assemble_LumpedSystem.c
ViewVC logotype

Contents of /trunk/finley/src/Assemble_LumpedSystem.c

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26