# Contents of /trunk-mpi-branch/finley/src/Assemble_PDE.c

Revision 1096 - (show annotations)
Mon Apr 16 22:59:33 2007 UTC (13 years, 1 month ago) by ksteube
File MIME type: text/plain
File size: 16825 byte(s)
```MPI implicit solver example run_simplesolve.py now compiling and
running successfully on one CPU of ess.

Adjusted SConscript, removed some debug print statements and removed
some partially implemented TRILINOS calls.

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

## Properties

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