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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 969 - (hide annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 1 month ago) by ksteube
Original Path: trunk/finley/src/Assemble_PDE.c
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 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 82
13 jgs 150
14 jgs 82 /**************************************************************/
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 jgs 150 /* Author: gross@access.edu.au */
43     /* Version: $Id$ */
44 jgs 82
45     /**************************************************************/
46    
47     #include "Assemble.h"
48     #include "Util.h"
49     #ifdef _OPENMP
50     #include <omp.h>
51     #endif
52    
53    
54     /**************************************************************/
55    
56 jgs 150 void Finley_Assemble_PDE(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
57 jgs 82 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
58    
59 gross 798 bool_t reducedIntegrationOrder=FALSE;
60 jgs 150 char error_msg[LenErrorMsg_MAX];
61 gross 798 Assemble_Parameters p;
62 jgs 82 double time0;
63 gross 798 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
64    
65 jgs 150 Finley_resetError();
66 jgs 82
67 ksteube 969 {
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 jgs 82 if (nodes==NULL || elements==NULL) return;
79     if (S==NULL && isEmpty(F)) return;
80    
81 gross 798 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 jgs 82
85 gross 798 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 jgs 82 }
88 gross 798
89     /* get the functionspace for this assemblage call */
90 jgs 123 type_t funcspace=UNKNOWN;
91 jgs 82 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 gross 798
100 jgs 82 /* check if all function spaces are the same */
101     if (! functionSpaceTypeEqual(funcspace,A) ) {
102 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A");
103 jgs 82 }
104     if (! functionSpaceTypeEqual(funcspace,B) ) {
105 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B");
106 jgs 82 }
107     if (! functionSpaceTypeEqual(funcspace,C) ) {
108 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C");
109 jgs 82 }
110     if (! functionSpaceTypeEqual(funcspace,D) ) {
111 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D");
112 jgs 82 }
113     if (! functionSpaceTypeEqual(funcspace,X) ) {
114 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X");
115 jgs 82 }
116     if (! functionSpaceTypeEqual(funcspace,Y) ) {
117 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y");
118 jgs 82 }
119 gross 798 if (! Finley_noError()) return;
120 jgs 82
121     /* check if all function spaces are the same */
122 gross 798 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 jgs 82
135 gross 798 /* 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 jgs 82 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
142 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
143 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
144 jgs 82 }
145    
146     if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
147 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
148 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
149 jgs 82 }
150    
151     if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
152 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
153 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
154 jgs 82 }
155    
156     if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
157 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
158 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
159 jgs 82 }
160    
161     if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
162 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
163 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
164 jgs 82 }
165    
166     if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
167 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
168 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
169 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
179 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
180 jgs 82 }
181     }
182     if (!isEmpty(B)) {
183     dimensions[0]=p.numDim;
184     if (!isDataPointShapeEqual(B,1,dimensions)) {
185 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
186 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
187 jgs 82 }
188     }
189     if (!isEmpty(C)) {
190     dimensions[0]=p.numDim;
191     if (!isDataPointShapeEqual(C,1,dimensions)) {
192 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
193 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
194 jgs 82 }
195     }
196     if (!isEmpty(D)) {
197     if (!isDataPointShapeEqual(D,0,dimensions)) {
198 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected.");
199 jgs 82 }
200     }
201     if (!isEmpty(X)) {
202     dimensions[0]=p.numDim;
203     if (!isDataPointShapeEqual(X,1,dimensions)) {
204 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
205 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
206 jgs 82 }
207     }
208     if (!isEmpty(Y)) {
209     if (!isDataPointShapeEqual(Y,0,dimensions)) {
210 gross 532 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected.");
211 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
221 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
222 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
230 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
231 jgs 82 }
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 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
239 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
240 jgs 82 }
241     }
242     if (!isEmpty(D)) {
243     dimensions[0]=p.numEqu;
244     dimensions[1]=p.numComp;
245     if (!isDataPointShapeEqual(D,2,dimensions)) {
246 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
247 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
248 jgs 82 }
249     }
250     if (!isEmpty(X)) {
251     dimensions[0]=p.numEqu;
252     dimensions[1]=p.numDim;
253     if (!isDataPointShapeEqual(X,2,dimensions)) {
254 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
255 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
256 jgs 82 }
257     }
258     if (!isEmpty(Y)) {
259     dimensions[0]=p.numEqu;
260     if (!isDataPointShapeEqual(Y,1,dimensions)) {
261 gross 532 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
262 jgs 150 Finley_setError(TYPE_ERROR,error_msg);
263 jgs 82 }
264     }
265     }
266 jgs 150 if (Finley_noError()) {
267 jgs 82 time0=Finley_timer();
268 gross 798 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 jgs 82 }
283 gross 798 } 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 jgs 82 }
355 jgs 149 #ifdef Finley_TRACE
356 jgs 82 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0);
357 jgs 149 #endif
358 jgs 82 }
359     }
360     /*
361     * $Log$
362 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 123 * 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 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
384 jgs 97 * *** empty log message ***
385 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
386     * some changes towards 64 integers in finley
387 jgs 82 *
388 jgs 123 * 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 jgs 97 *
391 jgs 82 *
392 jgs 123 *
393 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26