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

Revision 1353 - (hide annotations)
Thu Nov 22 06:05:35 2007 UTC (12 years, 3 months ago) by matt
File MIME type: text/plain
File size: 17109 byte(s)
```Reversed the removal of blocktimer code. This needs to be fixed properly for
the Windows stuff to work.

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

Properties

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