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

Revision 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (12 years, 10 months ago) by elspeth
Original Path: trunk/finley/src/Assemble_PDE.c
File MIME type: text/plain
File size: 16853 byte(s)
```Copyright added to more source files.

```
 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 char error_msg[LenErrorMsg_MAX]; 60 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Vol=NULL,*dvdV=NULL; 61 double time0; 62 dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q; 63 Assemble_Parameters p; 64 index_t *index_row=NULL,*index_col=NULL,color; 65 Finley_resetError(); 66 67 if (nodes==NULL || elements==NULL) return; 68 if (S==NULL && isEmpty(F)) return; 69 70 /* set all parameters in p*/ 71 Assemble_getAssembleParameters(nodes,elements,S,F,&p); 72 if (! Finley_noError()) return; 73 74 /* this function assumes NS=NN */ 75 if (p.NN!=p.NS) { 76 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: for Finley_Assemble_PDE numNodes and numShapes have to be identical."); 77 return; 78 } 79 if (p.numDim!=p.numElementDim) { 80 Finley_setError(SYSTEM_ERROR,"Finley_Assemble_PDE: Finley_Assemble_PDE accepts volume elements only."); 81 return; 82 } 83 /* get a functionspace */ 84 type_t funcspace=UNKNOWN; 85 updateFunctionSpaceType(funcspace,A); 86 updateFunctionSpaceType(funcspace,B); 87 updateFunctionSpaceType(funcspace,C); 88 updateFunctionSpaceType(funcspace,D); 89 updateFunctionSpaceType(funcspace,X); 90 updateFunctionSpaceType(funcspace,Y); 91 if (funcspace==UNKNOWN) return; /* all data are empty */ 92 93 /* check if all function spaces are the same */ 94 95 if (! functionSpaceTypeEqual(funcspace,A) ) { 96 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient A"); 97 } 98 if (! functionSpaceTypeEqual(funcspace,B) ) { 99 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient B"); 100 } 101 if (! functionSpaceTypeEqual(funcspace,C) ) { 102 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient C"); 103 } 104 if (! functionSpaceTypeEqual(funcspace,D) ) { 105 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient D"); 106 } 107 if (! functionSpaceTypeEqual(funcspace,X) ) { 108 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient X"); 109 } 110 if (! functionSpaceTypeEqual(funcspace,Y) ) { 111 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: unexpected function space type for coefficient Y"); 112 } 113 114 /* check if all function spaces are the same */ 115 116 if (! numSamplesEqual(A,p.numQuad,elements->numElements) ) { 117 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient A don't match (%d,%d)",p.numQuad,elements->numElements); 118 Finley_setError(TYPE_ERROR,error_msg); 119 } 120 121 if (! numSamplesEqual(B,p.numQuad,elements->numElements) ) { 122 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient B don't match (%d,%d)",p.numQuad,elements->numElements); 123 Finley_setError(TYPE_ERROR,error_msg); 124 } 125 126 if (! numSamplesEqual(C,p.numQuad,elements->numElements) ) { 127 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient C don't match (%d,%d)",p.numQuad,elements->numElements); 128 Finley_setError(TYPE_ERROR,error_msg); 129 } 130 131 if (! numSamplesEqual(D,p.numQuad,elements->numElements) ) { 132 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient D don't match (%d,%d)",p.numQuad,elements->numElements); 133 Finley_setError(TYPE_ERROR,error_msg); 134 } 135 136 if (! numSamplesEqual(X,p.numQuad,elements->numElements) ) { 137 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient X don't match (%d,%d)",p.numQuad,elements->numElements); 138 Finley_setError(TYPE_ERROR,error_msg); 139 } 140 141 if (! numSamplesEqual(Y,p.numQuad,elements->numElements) ) { 142 sprintf(error_msg,"Finley_Assemble_PDE: sample points of coefficient Y don't match (%d,%d)",p.numQuad,elements->numElements); 143 Finley_setError(TYPE_ERROR,error_msg); 144 } 145 146 /* check the dimensions: */ 147 148 if (p.numEqu==1 && p.numComp==1) { 149 if (!isEmpty(A)) { 150 dimensions[0]=p.numDim; 151 dimensions[1]=p.numDim; 152 if (!isDataPointShapeEqual(A,2,dimensions)) { 153 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A: illegal shape, expected shape (%d,%d)",dimensions[0],dimensions[1]); 154 Finley_setError(TYPE_ERROR,error_msg); 155 } 156 } 157 if (!isEmpty(B)) { 158 dimensions[0]=p.numDim; 159 if (!isDataPointShapeEqual(B,1,dimensions)) { 160 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B: illegal shape (%d,)",dimensions[0]); 161 Finley_setError(TYPE_ERROR,error_msg); 162 } 163 } 164 if (!isEmpty(C)) { 165 dimensions[0]=p.numDim; 166 if (!isDataPointShapeEqual(C,1,dimensions)) { 167 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,)",dimensions[0]); 168 Finley_setError(TYPE_ERROR,error_msg); 169 } 170 } 171 if (!isEmpty(D)) { 172 if (!isDataPointShapeEqual(D,0,dimensions)) { 173 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient D, rank 0 expected."); 174 } 175 } 176 if (!isEmpty(X)) { 177 dimensions[0]=p.numDim; 178 if (!isDataPointShapeEqual(X,1,dimensions)) { 179 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,",dimensions[0]); 180 Finley_setError(TYPE_ERROR,error_msg); 181 } 182 } 183 if (!isEmpty(Y)) { 184 if (!isDataPointShapeEqual(Y,0,dimensions)) { 185 Finley_setError(TYPE_ERROR,"Finley_Assemble_PDE: coefficient Y, rank 0 expected."); 186 } 187 } 188 } else { 189 if (!isEmpty(A)) { 190 dimensions[0]=p.numEqu; 191 dimensions[1]=p.numDim; 192 dimensions[2]=p.numComp; 193 dimensions[3]=p.numDim; 194 if (!isDataPointShapeEqual(A,4,dimensions)) { 195 sprintf(error_msg,"Finley_Assemble_PDE: coefficient A, expected shape (%d,%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2],dimensions[3]); 196 Finley_setError(TYPE_ERROR,error_msg); 197 } 198 } 199 if (!isEmpty(B)) { 200 dimensions[0]=p.numEqu; 201 dimensions[1]=p.numDim; 202 dimensions[2]=p.numComp; 203 if (!isDataPointShapeEqual(B,3,dimensions)) { 204 sprintf(error_msg,"Finley_Assemble_PDE: coefficient B, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]); 205 Finley_setError(TYPE_ERROR,error_msg); 206 } 207 } 208 if (!isEmpty(C)) { 209 dimensions[0]=p.numEqu; 210 dimensions[1]=p.numComp; 211 dimensions[2]=p.numDim; 212 if (!isDataPointShapeEqual(C,3,dimensions)) { 213 sprintf(error_msg,"Finley_Assemble_PDE: coefficient C, expected shape (%d,%d,%d)",dimensions[0],dimensions[1],dimensions[2]); 214 Finley_setError(TYPE_ERROR,error_msg); 215 } 216 } 217 if (!isEmpty(D)) { 218 dimensions[0]=p.numEqu; 219 dimensions[1]=p.numComp; 220 if (!isDataPointShapeEqual(D,2,dimensions)) { 221 sprintf(error_msg,"Finley_Assemble_PDE: coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]); 222 Finley_setError(TYPE_ERROR,error_msg); 223 } 224 } 225 if (!isEmpty(X)) { 226 dimensions[0]=p.numEqu; 227 dimensions[1]=p.numDim; 228 if (!isDataPointShapeEqual(X,2,dimensions)) { 229 sprintf(error_msg,"Finley_Assemble_PDE: coefficient X, expected shape (%d,%d)",dimensions[0],dimensions[1]); 230 Finley_setError(TYPE_ERROR,error_msg); 231 } 232 } 233 if (!isEmpty(Y)) { 234 dimensions[0]=p.numEqu; 235 if (!isDataPointShapeEqual(Y,1,dimensions)) { 236 sprintf(error_msg,"Finley_Assemble_PDE: coefficient Y, expected shape (%d,)",dimensions[0]); 237 Finley_setError(TYPE_ERROR,error_msg); 238 } 239 } 240 } 241 242 if (Finley_noError()) { 243 time0=Finley_timer(); 244 #pragma omp parallel private(index_col,index_row,EM_S,EM_F,V,dVdv,dSdV,Vol,dvdV,color,q) \ 245 firstprivate(elements,nodes,S,F,A,B,C,D,X,Y) 246 { 247 EM_S=EM_F=V=dVdv=dSdV=Vol=dvdV=NULL; 248 index_row=index_col=NULL; 249 250 /* allocate work arrays: */ 251 252 EM_S=(double*) THREAD_MEMALLOC(p.NN_row*p.NN_col*p.numEqu*p.numComp,double); 253 EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu,double); 254 V=(double*) THREAD_MEMALLOC(p.NN*p.numDim,double); 255 dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double); 256 dvdV=(double*) THREAD_MEMALLOC(p.numDim*p.numDim*p.numQuad,double); 257 dSdV=(double*) THREAD_MEMALLOC(p.NS_row*p.numQuad*p.numDim,double); 258 Vol=(double*) THREAD_MEMALLOC(p.numQuad,double); 259 index_col=(index_t*) THREAD_MEMALLOC(p.NN_col,index_t); 260 index_row=(index_t*) THREAD_MEMALLOC(p.NN_row,index_t); 261 262 if (! (Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(V) || Finley_checkPtr(index_col) || 263 Finley_checkPtr(index_row) || Finley_checkPtr(dVdv) || Finley_checkPtr(dSdV) || Finley_checkPtr(Vol) )) { 264 265 /* open loop over all colors: */ 266 for (color=elements->minColor;color<=elements->maxColor;color++) { 267 /* open loop over all elements: */ 268 #pragma omp for private(e) schedule(static) 269 for(e=0;enumElements;e++){ 270 if (elements->Color[e]==color) { 271 //============================ 272 for (q=0;qNodes[INDEX2(p.row_node[q],e,p.NN)]]; 273 /* gather V-coordinates of nodes into V: */ 274 Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V); 275 /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */ 276 Finley_Util_SmallMatMult(p.numDim,p.numDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv); 277 /* dvdV=invert(dVdv) inplace */ 278 Finley_Util_InvertSmallMat(p.numQuad,p.numDim,dVdv,dvdV,Vol); 279 /* calculate dSdV=DSDv*DvDV */ 280 Finley_Util_SmallMatSetMult(p.numQuad,p.NS_row,p.numDim,dSdV,p.numDim,p.referenceElement_row->dSdv,dvdV); 281 /* scale volume: */ 282 for (q=0;qQuadWeights[q]); 283 //============================ 284 285 /* integration for the stiffness matrix: */ 286 /* in order to optimze the number of operations the case of constants coefficience needs a bit more work */ 287 /* to make use of some symmetry. */ 288 289 if (S!=NULL) { 290 for (q=0;qS,dSdV,Vol,p.NN_row,EM_S, 294 getSampleData(A,e),isExpanded(A), 295 getSampleData(B,e),isExpanded(B), 296 getSampleData(C,e),isExpanded(C), 297 getSampleData(D,e),isExpanded(D)); 298 } else { 299 Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numDim,p.numQuad,p.numEqu,p.numComp, 300 p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_S, 301 getSampleData(A,e),isExpanded(A), 302 getSampleData(B,e),isExpanded(B), 303 getSampleData(C,e),isExpanded(C), 304 getSampleData(D,e),isExpanded(D)); 305 } 306 for (q=0;qNodes[INDEX2(p.col_node[q],e,p.NN)]]; 307 Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S); 308 } 309 if (!isEmpty(F)) { 310 for (q=0;qS,dSdV,Vol,p.NN_row,EM_F, 314 getSampleData(X,e),isExpanded(X), 315 getSampleData(Y,e),isExpanded(Y)); 316 } else { 317 Finley_Assemble_RHSMatrix_System(p.NS_row,p.numDim,p.numQuad, 318 p.numEqu,p.referenceElement_row->S,dSdV,Vol,p.NN_row,EM_F, 319 getSampleData(X,e),isExpanded(X), 320 getSampleData(Y,e),isExpanded(Y)); 321 } 322 /* add */ 323 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0)); 324 } 325 } 326 } 327 } 328 } 329 /* clean up */ 330 THREAD_MEMFREE(EM_S); 331 THREAD_MEMFREE(EM_F); 332 THREAD_MEMFREE(V); 333 THREAD_MEMFREE(dVdv); 334 THREAD_MEMFREE(dvdV); 335 THREAD_MEMFREE(dSdV); 336 THREAD_MEMFREE(Vol); 337 THREAD_MEMFREE(index_col); 338 THREAD_MEMFREE(index_row); 339 } 340 #ifdef Finley_TRACE 341 printf("timing: assemblage PDE: %.4e sec\n",Finley_timer()-time0); 342 #endif 343 } 344 } 345 /* 346 * \$Log\$ 347 * Revision 1.8 2005/09/15 03:44:21 jgs 348 * Merge of development branch dev-02 back to main trunk on 2005-09-15 349 * 350 * Revision 1.7 2005/09/01 03:31:35 jgs 351 * Merge of development branch dev-02 back to main trunk on 2005-09-01 352 * 353 * Revision 1.6 2005/08/12 01:45:42 jgs 354 * erge of development branch dev-02 back to main trunk on 2005-08-12 355 * 356 * Revision 1.5.2.3 2005/09/07 06:26:17 gross 357 * the solver from finley are put into the standalone package paso now 358 * 359 * Revision 1.5.2.2 2005/08/24 02:02:18 gross 360 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now. 361 * 362 * Revision 1.5.2.1 2005/08/03 08:54:27 gross 363 * contact element assemblage was called with wrong element table pointer 364 * 365 * Revision 1.5 2005/07/08 04:07:46 jgs 366 * Merge of development branch back to main trunk on 2005-07-08 367 * 368 * Revision 1.4 2004/12/15 07:08:32 jgs 369 * *** empty log message *** 370 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross 371 * some changes towards 64 integers in finley 372 * 373 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross 374 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now 375 * 376 * 377 * 378 */

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26