# Contents of /trunk/finley/src/Assemble_PDE.c

Revision 2156 - (show annotations)
Mon Dec 15 05:09:02 2008 UTC (10 years, 11 months ago) by gross
File MIME type: text/plain
File size: 17076 byte(s)
```some modifications to the iterative solver to make them easier to use.
There are also improved versions of the Darcy flux solver and the incompressible solver.

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

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26