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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * 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
14
15 /**************************************************************/
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 #include "esysUtils/blocktimer.h"
46 #ifdef _OPENMP
47 #include <omp.h>
48 #endif
49
50
51 /**************************************************************/
52
53 void Dudley_Assemble_PDE(Dudley_NodeFile* nodes,Dudley_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 double blocktimer_start = blocktimer_time();
63
64 Dudley_resetError();
65
66 {
67 #ifdef Paso_MPI
68 int iam, numCPUs;
69 iam = elements->elementDistribution->MPIInfo->rank;
70 numCPUs = elements->elementDistribution->MPIInfo->size;
71 #endif
72 }
73
74 if (nodes==NULL || elements==NULL) return;
75 if (S==NULL && isEmpty(F)) return;
76
77 if (isEmpty(F) && !isEmpty(X) && !isEmpty(F)) {
78 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: right hand side coefficients are non-zero bat no right hand side vector given.");
79 }
80
81 if (S==NULL && !isEmpty(A) && !isEmpty(B) && !isEmpty(C) && !isEmpty(D)) {
82 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: coefficients are non-zero but no matrix is given.");
83 }
84
85 /* get the functionspace for this assemblage call */
86 funcspace=UNKNOWN;
87 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 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient A");
98 }
99 if (! functionSpaceTypeEqual(funcspace,B) ) {
100 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient B");
101 }
102 if (! functionSpaceTypeEqual(funcspace,C) ) {
103 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient C");
104 }
105 if (! functionSpaceTypeEqual(funcspace,D) ) {
106 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient D");
107 }
108 if (! functionSpaceTypeEqual(funcspace,X) ) {
109 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient X");
110 }
111 if (! functionSpaceTypeEqual(funcspace,Y) ) {
112 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: unexpected function space type for coefficient Y");
113 }
114 if (! Dudley_noError()) return;
115
116 /* check if all function spaces are the same */
117 if (funcspace==DUDLEY_ELEMENTS) {
118 reducedIntegrationOrder=FALSE;
119 } else if (funcspace==DUDLEY_FACE_ELEMENTS) {
120 reducedIntegrationOrder=FALSE;
121 } else if (funcspace==DUDLEY_REDUCED_ELEMENTS) {
122 reducedIntegrationOrder=TRUE;
123 } else if (funcspace==DUDLEY_REDUCED_FACE_ELEMENTS) {
124 reducedIntegrationOrder=TRUE;
125 } else {
126 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: assemblage failed because of illegal function space.");
127 }
128 if (! Dudley_noError()) return;
129
130 /* set all parameters in p*/
131 Assemble_getAssembleParameters(nodes,elements,S,F, reducedIntegrationOrder, &p);
132 if (! Dudley_noError()) return;
133
134 /* check if all function spaces are the same */
135
136 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 Dudley_setError(TYPE_ERROR,error_msg);
139 }
140
141 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 Dudley_setError(TYPE_ERROR,error_msg);
144 }
145
146 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 Dudley_setError(TYPE_ERROR,error_msg);
149 }
150
151 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 Dudley_setError(TYPE_ERROR,error_msg);
154 }
155
156 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 Dudley_setError(TYPE_ERROR,error_msg);
159 }
160
161 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 Dudley_setError(TYPE_ERROR,error_msg);
164 }
165
166 /* check the dimensions: */
167
168 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 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 }
176 }
177 if (!isEmpty(B)) {
178 dimensions[0]=p.numDim;
179 if (!isDataPointShapeEqual(B,1,dimensions)) {
180 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]);
181 Dudley_setError(TYPE_ERROR,error_msg);
182 }
183 }
184 if (!isEmpty(C)) {
185 dimensions[0]=p.numDim;
186 if (!isDataPointShapeEqual(C,1,dimensions)) {
187 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]);
188 Dudley_setError(TYPE_ERROR,error_msg);
189 }
190 }
191 if (!isEmpty(D)) {
192 if (!isDataPointShapeEqual(D,0,dimensions)) {
193 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: coefficient D, rank 0 expected.");
194 }
195 }
196 if (!isEmpty(X)) {
197 dimensions[0]=p.numDim;
198 if (!isDataPointShapeEqual(X,1,dimensions)) {
199 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]);
200 Dudley_setError(TYPE_ERROR,error_msg);
201 }
202 }
203 if (!isEmpty(Y)) {
204 if (!isDataPointShapeEqual(Y,0,dimensions)) {
205 Dudley_setError(TYPE_ERROR,"Dudley_Assemble_PDE: coefficient Y, rank 0 expected.");
206 }
207 }
208 } 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 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 }
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 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 }
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 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 }
236 }
237 if (!isEmpty(D)) {
238 dimensions[0]=p.numEqu;
239 dimensions[1]=p.numComp;
240 if (!isDataPointShapeEqual(D,2,dimensions)) {
241 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
242 Dudley_setError(TYPE_ERROR,error_msg);
243 }
244 }
245 if (!isEmpty(X)) {
246 dimensions[0]=p.numEqu;
247 dimensions[1]=p.numDim;
248 if (!isDataPointShapeEqual(X,2,dimensions)) {
249 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]);
250 Dudley_setError(TYPE_ERROR,error_msg);
251 }
252 }
253 if (!isEmpty(Y)) {
254 dimensions[0]=p.numEqu;
255 if (!isDataPointShapeEqual(Y,1,dimensions)) {
256 sprintf(error_msg,"Dudley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]);
257 Dudley_setError(TYPE_ERROR,error_msg);
258 }
259 }
260 }
261 if (Dudley_noError()) {
262 time0=Dudley_timer();
263 if (p.numEqu == p. numComp) {
264 if (p.numEqu > 1) {
265 /* system of PDESs */
266 if (p.numDim==3) {
267 Dudley_Assemble_PDE_System2_3D(p,elements,S,F,A,B,C,D,X,Y);
268 } else if (p.numDim==2) {
269 Dudley_Assemble_PDE_System2_2D(p,elements,S,F,A,B,C,D,X,Y);
270 } else {
271 Dudley_setError(VALUE_ERROR,"Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
272 }
273 } else {
274 /* single PDES */
275 if (p.numDim==3) {
276 Dudley_Assemble_PDE_Single2_3D(p,elements,S,F,A,B,C,D,X,Y);
277 } else if (p.numDim==2) {
278 Dudley_Assemble_PDE_Single2_2D(p,elements,S,F,A,B,C,D,X,Y);
279 } else {
280 Dudley_setError(VALUE_ERROR,"Dudley_Assemble_PDE supports spatial dimensions 2 and 3 only.");
281 }
282 }
283 } else {
284 Dudley_setError(VALUE_ERROR,"Dudley_Assemble_PDE requires number of equations == number of solutions .");
285 }
286 }
287 blocktimer_increment("Dudley_Assemble_PDE()", blocktimer_start);
288 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26