/[escript]/branches/domexper/dudley/src/Assemble_PDE.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Assemble_PDE.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3205 - (hide annotations)
Fri Sep 24 00:30:43 2010 UTC (9 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 11185 byte(s)
Removing -total from field names

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
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 phornby 2078 #include "esysUtils/blocktimer.h"
46 jgs 82 #ifdef _OPENMP
47     #include <omp.h>
48     #endif
49    
50    
51     /**************************************************************/
52    
53 jfenwick 3086 void Dudley_Assemble_PDE(Dudley_NodeFile* nodes,Dudley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F,
54 jgs 82 escriptDataC* A, escriptDataC* B, escriptDataC* C, escriptDataC* D, escriptDataC* X, escriptDataC* Y ) {
55    
56 gross 798 bool_t reducedIntegrationOrder=FALSE;
57 jgs 150 char error_msg[LenErrorMsg_MAX];
58 gross 798 Assemble_Parameters p;
59 jgs 82 double time0;
60 gross 798 dim_t dimensions[ESCRIPT_MAX_DATA_RANK];
61 gross 1028 type_t funcspace;
62 matt 1353 double blocktimer_start = blocktimer_time();
63 gross 798
64 jfenwick 3086 Dudley_resetError();
65 jgs 82
66 ksteube 1312 {
67 phornby 1628 #ifdef Paso_MPI
68 ksteube 1312 int iam, numCPUs;
69     iam = elements->elementDistribution->MPIInfo->rank;
70     numCPUs = elements->elementDistribution->MPIInfo->size;
71     #endif
72     }
73    
74 jgs 82 if (nodes==NULL || elements==NULL) return;
75     if (S==NULL && isEmpty(F)) return;
76    
77 matt 1352 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
78 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
79 gross 798 }
80 jgs 82
81 matt 1352 if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
82 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
83 jgs 82 }
84 gross 798
85     /* get the functionspace for this assemblage call */
86 gross 1028 funcspace=UNKNOWN;
87 jgs 82 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 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient A");
98 jgs 82 }
99     if (! functionSpaceTypeEqual(funcspace,B) ) {
100 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient B");
101 jgs 82 }
102     if (! functionSpaceTypeEqual(funcspace,C) ) {
103 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient C");
104 jgs 82 }
105     if (! functionSpaceTypeEqual(funcspace,D) ) {
106 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient D");
107 jgs 82 }
108     if (! functionSpaceTypeEqual(funcspace,X) ) {
109 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient X");
110 jgs 82 }
111     if (! functionSpaceTypeEqual(funcspace,Y) ) {
112 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient Y");
113 jgs 82 }
114 jfenwick 3086 if (! Dudley_noError()) return;
115 jgs 82
116     /* check if all function spaces are the same */
117 jfenwick 3086 if (funcspace==DUDLEY_ELEMENTS) {
118 gross 798 reducedIntegrationOrder=FALSE;
119 jfenwick 3086 } else if (funcspace==DUDLEY_FACE_ELEMENTS) {
120 gross 798 reducedIntegrationOrder=FALSE;
121 jfenwick 3086 } else if (funcspace==DUDLEY_REDUCED_ELEMENTS) {
122 gross 1072 reducedIntegrationOrder=TRUE;
123 jfenwick 3086 } else if (funcspace==DUDLEY_REDUCED_FACE_ELEMENTS) {
124 gross 1072 reducedIntegrationOrder=TRUE;
125 gross 798 } else {
126 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: assemblage failed because of illegal function space.");
127 gross 798 }
128 jfenwick 3086 if (! Dudley_noError()) return;
129 jgs 82
130 gross 798 /* set all parameters in p*/
131     Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
132 jfenwick 3086 if (! Dudley_noError()) return;
133 gross 798
134     /* check if all function spaces are the same */
135    
136 jfenwick 3205 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) {
137     sprintf(error_msg,"Dudley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements);
138 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
139 jgs 82 }
140    
141 jfenwick 3205 if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) {
142     sprintf(error_msg,"Dudley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements);
143 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
144 jgs 82 }
145    
146 jfenwick 3205 if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) {
147     sprintf(error_msg,"Dudley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements);
148 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
149 jgs 82 }
150    
151 jfenwick 3205 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) {
152     sprintf(error_msg,"Dudley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements);
153 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
154 jgs 82 }
155    
156 jfenwick 3205 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) {
157     sprintf(error_msg,"Dudley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements);
158 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
159 jgs 82 }
160    
161 jfenwick 3205 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) {
162     sprintf(error_msg,"Dudley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements);
163 jfenwick 3086 Dudley_setError(TYPE_ERROR,error_msg);
164 jgs 82 }
165    
166     /* check the dimensions: */
167 matt 1352
168 jgs 82 if (p.numEqu==1 && p.numComp==1) {
169     if (!isEmpty(A)) {
170     dimensions[0]=p.numDim;
171     dimensions[1]=p.numDim;
172     if (!isDataPointShapeEqual(A,2,dimensions)) {
173 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]);
174     Dudley_setError(TYPE_ERROR,error_msg);
175 jgs 82 }
176     }
177     if (!isEmpty(B)) {
178     dimensions[0]=p.numDim;
179     if (!isDataPointShapeEqual(B,1,dimensions)) {
180 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
181     Dudley_setError(TYPE_ERROR,error_msg);
182 jgs 82 }
183     }
184     if (!isEmpty(C)) {
185     dimensions[0]=p.numDim;
186     if (!isDataPointShapeEqual(C,1,dimensions)) {
187 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
188     Dudley_setError(TYPE_ERROR,error_msg);
189 jgs 82 }
190     }
191     if (!isEmpty(D)) {
192     if (!isDataPointShapeEqual(D,0,dimensions)) {
193 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: coefficient D, rank 0 expected.");
194 jgs 82 }
195     }
196     if (!isEmpty(X)) {
197     dimensions[0]=p.numDim;
198     if (!isDataPointShapeEqual(X,1,dimensions)) {
199 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
200     Dudley_setError(TYPE_ERROR,error_msg);
201 jgs 82 }
202     }
203     if (!isEmpty(Y)) {
204     if (!isDataPointShapeEqual(Y,0,dimensions)) {
205 jfenwick 3086 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: coefficient Y, rank 0 expected.");
206 jgs 82 }
207 matt 1352 }
208 jgs 82 } else {
209     if (!isEmpty(A)) {
210     dimensions[0]=p.numEqu;
211     dimensions[1]=p.numDim;
212     dimensions[2]=p.numComp;
213     dimensions[3]=p.numDim;
214     if (!isDataPointShapeEqual(A,4,dimensions)) {
215 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]);
216     Dudley_setError(TYPE_ERROR,error_msg);
217 jgs 82 }
218     }
219     if (!isEmpty(B)) {
220     dimensions[0]=p.numEqu;
221     dimensions[1]=p.numDim;
222     dimensions[2]=p.numComp;
223     if (!isDataPointShapeEqual(B,3,dimensions)) {
224 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
225     Dudley_setError(TYPE_ERROR,error_msg);
226 jgs 82 }
227     }
228     if (!isEmpty(C)) {
229     dimensions[0]=p.numEqu;
230     dimensions[1]=p.numComp;
231     dimensions[2]=p.numDim;
232     if (!isDataPointShapeEqual(C,3,dimensions)) {
233 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]);
234     Dudley_setError(TYPE_ERROR,error_msg);
235 jgs 82 }
236     }
237     if (!isEmpty(D)) {
238     dimensions[0]=p.numEqu;
239     dimensions[1]=p.numComp;
240     if (!isDataPointShapeEqual(D,2,dimensions)) {
241 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
242     Dudley_setError(TYPE_ERROR,error_msg);
243 jgs 82 }
244     }
245     if (!isEmpty(X)) {
246     dimensions[0]=p.numEqu;
247     dimensions[1]=p.numDim;
248     if (!isDataPointShapeEqual(X,2,dimensions)) {
249 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
250     Dudley_setError(TYPE_ERROR,error_msg);
251 jgs 82 }
252     }
253     if (!isEmpty(Y)) {
254     dimensions[0]=p.numEqu;
255     if (!isDataPointShapeEqual(Y,1,dimensions)) {
256 jfenwick 3086 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
257     Dudley_setError(TYPE_ERROR,error_msg);
258 jgs 82 }
259     }
260     }
261 jfenwick 3086 if (Dudley_noError()) {
262     time0=Dudley_timer();
263 gross 798 if (p.numEqu == p. numComp) {
264     if (p.numEqu > 1) {
265     /* system of PDESs */
266     if (p.numDim==3) {
267 jfenwick 3086 Dudley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
268 gross 798 } else if (p.numDim==2) {
269 jfenwick 3086 Dudley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
270 gross 798 } else {
271 jfenwick 3156 Dudley_setError(VALUE_ERROR,"Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
272 gross 798 }
273     } else {
274     /* single PDES */
275     if (p.numDim==3) {
276 jfenwick 3086 Dudley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
277 gross 798 } else if (p.numDim==2) {
278 jfenwick 3086 Dudley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
279 gross 798 } else {
280 jfenwick 3156 Dudley_setError(VALUE_ERROR,"Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
281 gross 798 }
282     }
283     } else {
284 jfenwick 3086 Dudley_setError(VALUE_ERROR,"Dudley_Assemble_PDE requires number of equations == number of solutions .");
285 jgs 82 }
286     }
287 jfenwick 3086 blocktimer_increment("Dudley_Assemble_PDE()", blocktimer_start);
288 jgs 82 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26