![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 17095
- Timestamp:
- 01/10/14 13:42:13 (11 years ago)
- Location:
- issm/trunk-jpl/src/c/classes/Elements
- Files:
-
- 4 edited
-
PentaRef.cpp (modified) (1 diff)
-
PentaRef.h (modified) (1 diff)
-
TriaRef.cpp (modified) (2 diffs)
-
TriaRef.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
r17070 r17095 56 56 57 57 /*Reference Element numerics*/ 58 /*FUNCTION PentaRef::GetBSSAHO {{{*/59 void PentaRef::GetBSSAHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){60 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.61 * For node i, Bi can be expressed in the actual coordinate system62 * by:63 * Bi=[ dh/dx 0 ]64 * [ 0 dh/dy ]65 * [ 1/2*dh/dy 1/2*dh/dx ]66 * where h is the interpolation function for node i.67 *68 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)69 */70 71 IssmDouble dbasis[3][NUMNODESP1];72 73 /*Get dbasis in actual coordinate system: */74 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);75 76 /*Build B: */77 for(int i=0;i<NUMNODESP1;i++){78 B[NDOF2*NUMNODESP1*0+NDOF2*i+0] = dbasis[0][i];79 B[NDOF2*NUMNODESP1*0+NDOF2*i+1] = 0.;80 81 B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;82 B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];83 84 B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];85 B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];86 }87 }88 /*}}}*/89 /*FUNCTION PentaRef::GetBSSAFS{{{*/90 void PentaRef::GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){91 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.92 * For node i, Bi can be expressed in the actual coordinate system93 * by:94 * Bi=[ dh/dx 0 0 0 ]95 * [ 0 dh/dy 0 0 ]96 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ]97 * [ 0 0 0 h ]98 * where h is the interpolation function for node i.99 *100 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)101 */102 103 int i;104 IssmDouble dbasismini[3][NUMNODESP1b];105 IssmDouble basis[NUMNODESP1];106 107 /*Get dbasis in actual coordinate system: */108 GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);109 GetNodalFunctionsP1(basis,gauss);110 111 /*Build B: */112 for(i=0;i<NUMNODESP1;i++){113 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+0] = dbasismini[0][i];114 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+1] = 0.;115 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+2] = 0.;116 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+0] = 0.;117 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+1] = dbasismini[1][i];118 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+2] = 0.;119 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+0] = 0.5*dbasismini[1][i];120 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+1] = 0.5*dbasismini[0][i];121 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+2] = 0.;122 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+0] = 0.;123 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+1] = 0.;124 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+2] = 0.;125 }126 for(i=0;i<1;i++){127 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+0] = 0.;128 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+1] = 0.;129 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+2] = 0.;130 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+0] = 0.;131 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+1] = 0.;132 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+2] = 0.;133 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+0] = 0.;134 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+1] = 0.;135 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+2] = 0.;136 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+0] = 0.;137 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+1] = 0.;138 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+2] = 0.;139 }140 141 for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function142 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NUMNODESP1b*NDOF3+i] = 0;143 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NUMNODESP1b*NDOF3+i] = 0;144 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NUMNODESP1b*NDOF3+i] = 0;145 B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NUMNODESP1b*NDOF3+i] = basis[i];146 }147 }148 /*}}}*/149 /*FUNCTION PentaRef::GetBHO {{{*/150 void PentaRef::GetBHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){151 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.152 * For node i, Bi can be expressed in the actual coordinate system153 * by:154 * Bi=[ dh/dx 0 ]155 * [ 0 dh/dy ]156 * [ 1/2*dh/dy 1/2*dh/dx ]157 * [ 1/2*dh/dz 0 ]158 * [ 0 1/2*dh/dz ]159 * where h is the interpolation function for node i.160 *161 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)162 */163 164 /*Fetch number of nodes for this finite element*/165 int numnodes = this->NumberofNodes();166 167 /*Get nodal functions derivatives*/168 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);169 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);170 171 /*Build B: */172 for(int i=0;i<numnodes;i++){173 B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];174 B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;175 B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;176 B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];177 B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];178 B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];179 B[NDOF2*numnodes*3+NDOF2*i+0] = .5*dbasis[2*numnodes+i];180 B[NDOF2*numnodes*3+NDOF2*i+1] = 0.;181 B[NDOF2*numnodes*4+NDOF2*i+0] = 0.;182 B[NDOF2*numnodes*4+NDOF2*i+1] = .5*dbasis[2*numnodes+i];183 }184 185 /*Clean-up*/186 xDelete<IssmDouble>(dbasis);187 }188 /*}}}*/189 /*FUNCTION PentaRef::GetBprimeHO {{{*/190 void PentaRef::GetBprimeHO(IssmDouble* B,IssmDouble* xyz_list,Gauss* gauss){191 /*Compute B prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.192 * For node i, Bi can be expressed in the actual coordinate system193 * by:194 * Bi=[ 2*dh/dx dh/dy ]195 * [ dh/dx 2*dh/dy ]196 * [ dh/dy dh/dx ]197 * [ dh/dz 0 ]198 * [ 0 dh/dz ]199 * where h is the interpolation function for node i.200 *201 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)202 */203 204 /*Fetch number of nodes for this finite element*/205 int numnodes = this->NumberofNodes();206 207 /*Get nodal functions derivatives*/208 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);209 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);210 211 /*Build BPrime: */212 for(int i=0;i<numnodes;i++){213 B[NDOF2*numnodes*0+NDOF2*i+0]=2.*dbasis[0*numnodes+i];214 B[NDOF2*numnodes*0+NDOF2*i+1]=dbasis[1*numnodes+i];215 B[NDOF2*numnodes*1+NDOF2*i+0]=dbasis[0*numnodes+i];216 B[NDOF2*numnodes*1+NDOF2*i+1]=2.*dbasis[1*numnodes+i];217 B[NDOF2*numnodes*2+NDOF2*i+0]=dbasis[1*numnodes+i];218 B[NDOF2*numnodes*2+NDOF2*i+1]=dbasis[0*numnodes+i];219 B[NDOF2*numnodes*3+NDOF2*i+0]=dbasis[2*numnodes+i];220 B[NDOF2*numnodes*3+NDOF2*i+1]=0.;221 B[NDOF2*numnodes*4+NDOF2*i+0]=0.;222 B[NDOF2*numnodes*4+NDOF2*i+1]=dbasis[2*numnodes+i];223 }224 225 /*Clean-up*/226 xDelete<IssmDouble>(dbasis);227 }228 /*}}}*/229 /*FUNCTION PentaRef::GetBprimeSSAFS{{{*/230 void PentaRef::GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){231 /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.232 * For node i, Bprimei can be expressed in the actual coordinate system233 * by:234 * Bprimei=[ 2*dh/dx dh/dy 0 0 ]235 * [ dh/dx 2*dh/dy 0 0 ]236 * [ dh/dy dh/dx 0 0 ]237 * where h is the interpolation function for node i.238 *239 * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)240 */241 242 int i;243 IssmDouble dbasismini[3][NUMNODESP1b];244 245 /*Get dbasis in actual coordinate system: */246 GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);247 248 /*Build Bprime: */249 for(i=0;i<NUMNODESP1;i++){250 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+0] = 2.*dbasismini[0][i];251 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+1] = dbasismini[1][i];252 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+2] = 0.;253 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+0] = dbasismini[0][i];254 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+1] = 2.*dbasismini[1][i];255 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+2] = 0.;256 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+0] = dbasismini[1][i];257 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+1] = dbasismini[0][i];258 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+2] = 0.;259 }260 261 for(i=0;i<1;i++){ //Add zeros for the bubble function262 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+0] = 0.;263 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+1] = 0.;264 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+2] = 0.;265 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+0] = 0.;266 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+1] = 0.;267 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+2] = 0.;268 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+0] = 0.;269 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+1] = 0.;270 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+2] = 0.;271 }272 273 for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function274 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NUMNODESP1b*NDOF3+i] = 0.;275 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NUMNODESP1b*NDOF3+i] = 0.;276 Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NUMNODESP1b*NDOF3+i] = 0.;277 }278 279 }280 /*}}}*/281 /*FUNCTION PentaRef::GetBFSstrainrate {{{*/282 void PentaRef::GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){283 284 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.285 * For node i, Bi can be expressed in the actual coordinate system286 * by: Bi=[ dh/dx 0 0 ]287 * [ 0 dh/dy 0 ]288 * [ 0 0 dh/dy ]289 * [ 1/2*dh/dy 1/2*dh/dx 0 ]290 * [ 1/2*dh/dz 0 1/2*dh/dx ]291 * [ 0 1/2*dh/dz 1/2*dh/dy ]292 * where h is the interpolation function for node i.293 * Same thing for Bb except the last column that does not exist.294 */295 296 /*Fetch number of nodes for this finite element*/297 int numnodes = this->NumberofNodes();298 299 /*Get nodal functions derivatives*/300 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);301 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);302 303 /*Build B: */304 for(int i=0;i<numnodes;i++){305 B[3*numnodes*0+3*i+0] = dbasis[0*numnodes+i+0];306 B[3*numnodes*0+3*i+1] = 0.;307 B[3*numnodes*0+3*i+2] = 0.;308 309 B[3*numnodes*1+3*i+0] = 0.;310 B[3*numnodes*1+3*i+1] = dbasis[1*numnodes+i+0];311 B[3*numnodes*1+3*i+2] = 0.;312 313 B[3*numnodes*2+3*i+0] = 0.;314 B[3*numnodes*2+3*i+1] = 0.;315 B[3*numnodes*2+3*i+2] = dbasis[2*numnodes+i+0];316 317 B[3*numnodes*3+3*i+0] = .5*dbasis[1*numnodes+i+0];318 B[3*numnodes*3+3*i+1] = .5*dbasis[0*numnodes+i+0];319 B[3*numnodes*3+3*i+2] = 0.;320 321 B[3*numnodes*4+3*i+0] = .5*dbasis[2*numnodes+i+0];322 B[3*numnodes*4+3*i+1] = 0.;323 B[3*numnodes*4+3*i+2] = .5*dbasis[0*numnodes+i+0];324 325 B[3*numnodes*5+3*i+0] = 0.;326 B[3*numnodes*5+3*i+1] = .5*dbasis[2*numnodes+i+0];327 B[3*numnodes*5+3*i+2] = .5*dbasis[1*numnodes+i+0];328 }329 330 /*Clean up*/331 xDelete<IssmDouble>(dbasis);332 }333 /*}}}*/334 /*FUNCTION PentaRef::GetBFS {{{*/335 void PentaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){336 337 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.338 * For node i, Bvi can be expressed in the actual coordinate system339 * by: Bvi=[ dh/dx 0 0 ]340 * [ 0 dh/dy 0 ]341 * [ 0 0 dh/dz ]342 * [ 1/2*dh/dy 1/2*dh/dx 0 ]343 * [ 1/2*dh/dz 0 1/2*dh/dx ]344 * [ 0 1/2*dh/dz 1/2*dh/dy ]345 * [ 0 0 0 ]346 * [ dh/dx dh/dy dh/dz ]347 *348 * by: Bpi=[ 0 ]349 * [ 0 ]350 * [ 0 ]351 * [ 0 ]352 * [ 0 ]353 * [ 0 ]354 * [ h ]355 * [ 0 ]356 * where h is the interpolation function for node i.357 * Same thing for Bb except the last column that does not exist.358 */359 360 /*Fetch number of nodes for this finite element*/361 int pnumnodes = this->NumberofNodesPressure();362 int vnumnodes = this->NumberofNodesVelocity();363 364 /*Get nodal functions derivatives*/365 IssmDouble* vdbasis=xNew<IssmDouble>(3*vnumnodes);366 IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);367 GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);368 GetNodalFunctionsPressure(pbasis,gauss);369 370 /*Build B: */371 for(int i=0;i<vnumnodes;i++){372 B[(3*vnumnodes+pnumnodes)*0+3*i+0] = vdbasis[0*vnumnodes+i];373 B[(3*vnumnodes+pnumnodes)*0+3*i+1] = 0.;374 B[(3*vnumnodes+pnumnodes)*0+3*i+2] = 0.;375 B[(3*vnumnodes+pnumnodes)*1+3*i+0] = 0.;376 B[(3*vnumnodes+pnumnodes)*1+3*i+1] = vdbasis[1*vnumnodes+i];377 B[(3*vnumnodes+pnumnodes)*1+3*i+2] = 0.;378 B[(3*vnumnodes+pnumnodes)*2+3*i+0] = 0.;379 B[(3*vnumnodes+pnumnodes)*2+3*i+1] = 0.;380 B[(3*vnumnodes+pnumnodes)*2+3*i+2] = vdbasis[2*vnumnodes+i];381 B[(3*vnumnodes+pnumnodes)*3+3*i+0] = .5*vdbasis[1*vnumnodes+i];382 B[(3*vnumnodes+pnumnodes)*3+3*i+1] = .5*vdbasis[0*vnumnodes+i];383 B[(3*vnumnodes+pnumnodes)*3+3*i+2] = 0.;384 B[(3*vnumnodes+pnumnodes)*4+3*i+0] = .5*vdbasis[2*vnumnodes+i];385 B[(3*vnumnodes+pnumnodes)*4+3*i+1] = 0.;386 B[(3*vnumnodes+pnumnodes)*4+3*i+2] = .5*vdbasis[0*vnumnodes+i];387 B[(3*vnumnodes+pnumnodes)*5+3*i+0] = 0.;388 B[(3*vnumnodes+pnumnodes)*5+3*i+1] = .5*vdbasis[2*vnumnodes+i];389 B[(3*vnumnodes+pnumnodes)*5+3*i+2] = .5*vdbasis[1*vnumnodes+i];390 B[(3*vnumnodes+pnumnodes)*6+3*i+0] = 0.;391 B[(3*vnumnodes+pnumnodes)*6+3*i+1] = 0.;392 B[(3*vnumnodes+pnumnodes)*6+3*i+2] = 0.;393 B[(3*vnumnodes+pnumnodes)*7+3*i+0] = vdbasis[0*vnumnodes+i];394 B[(3*vnumnodes+pnumnodes)*7+3*i+1] = vdbasis[1*vnumnodes+i];395 B[(3*vnumnodes+pnumnodes)*7+3*i+2] = vdbasis[2*vnumnodes+i];396 }397 for(int i=0;i<pnumnodes;i++){398 B[(3*vnumnodes+pnumnodes)*0+(3*vnumnodes)+i] = 0.;399 B[(3*vnumnodes+pnumnodes)*1+(3*vnumnodes)+i] = 0.;400 B[(3*vnumnodes+pnumnodes)*2+(3*vnumnodes)+i] = 0.;401 B[(3*vnumnodes+pnumnodes)*3+(3*vnumnodes)+i] = 0.;402 B[(3*vnumnodes+pnumnodes)*4+(3*vnumnodes)+i] = 0.;403 B[(3*vnumnodes+pnumnodes)*5+(3*vnumnodes)+i] = 0.;404 B[(3*vnumnodes+pnumnodes)*6+(3*vnumnodes)+i] = pbasis[i];405 B[(3*vnumnodes+pnumnodes)*7+(3*vnumnodes)+i] = 0.;406 }407 408 /*Clean up*/409 xDelete<IssmDouble>(vdbasis);410 xDelete<IssmDouble>(pbasis);411 }412 /*}}}*/413 /*FUNCTION PentaRef::GetBFSGLS {{{*/414 void PentaRef::GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){415 416 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.417 * For node i, Bi can be expressed in the actual coordinate system418 * by: Bi=[ dh/dx 0 0 0 ]419 * [ 0 dh/dy 0 0 ]420 * [ 0 0 dh/dy 0 ]421 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ]422 * [ 1/2*dh/dz 0 1/2*dh/dx 0 ]423 * [ 0 1/2*dh/dz 1/2*dh/dy 0 ]424 * [ 0 0 0 h ]425 * [ dh/dx dh/dy dh/dz 0 ]426 * where h is the interpolation function for node i.427 * Same thing for Bb except the last column that does not exist.428 */429 430 int i;431 IssmDouble dbasis[3][NUMNODESP1];432 IssmDouble basis[NUMNODESP1];433 434 /*Get dbasismini in actual coordinate system: */435 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);436 GetNodalFunctionsP1(&basis[0], gauss);437 438 /*Build B: */439 for(i=0;i<NUMNODESP1;i++){440 B[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];441 B[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;442 B[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;443 B[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;444 B[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];445 B[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;446 B[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;447 B[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;448 B[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];449 B[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = .5*dbasis[1][i];450 B[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = .5*dbasis[0][i];451 B[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;452 B[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = .5*dbasis[2][i];453 B[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;454 B[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = .5*dbasis[0][i];455 B[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;456 B[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = .5*dbasis[2][i];457 B[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = .5*dbasis[1][i];458 B[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = 0.;459 B[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = 0.;460 B[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = 0.;461 B[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = dbasis[0][i];462 B[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = dbasis[1][i];463 B[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = dbasis[2][i];464 }465 466 for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function467 B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;468 B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;469 B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;470 B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;471 B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;472 B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;473 B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = basis[i];474 B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = 0.;475 }476 477 }478 /*}}}*/479 /*FUNCTION PentaRef::GetBprimeFS {{{*/480 void PentaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss){481 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.482 * For node i, Bi' can be expressed in the actual coordinate system483 * by:484 * Bvi' = [ dh/dx 0 0 ]485 * [ 0 dh/dy 0 ]486 * [ 0 0 dh/dz ]487 * [ dh/dy dh/dx 0 ]488 * [ dh/dz 0 dh/dx ]489 * [ 0 dh/dz dh/dy ]490 * [ dh/dx dh/dy dh/dz ]491 * [ 0 0 0 ]492 *493 * by: Bpi=[ 0 ]494 * [ 0 ]495 * [ 0 ]496 * [ 0 ]497 * [ 0 ]498 * [ 0 ]499 * [ 0 ]500 * [ h ]501 * where h is the interpolation function for node i.502 */503 504 /*Fetch number of nodes for this finite element*/505 int pnumnodes = this->NumberofNodesPressure();506 int vnumnodes = this->NumberofNodesVelocity();507 508 /*Get nodal functions derivatives*/509 IssmDouble* vdbasis=xNew<IssmDouble>(3*vnumnodes);510 IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);511 GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);512 GetNodalFunctionsPressure(pbasis,gauss);513 514 /*Build B_prime: */515 for(int i=0;i<vnumnodes;i++){516 B_prime[(3*vnumnodes+pnumnodes)*0+3*i+0] = vdbasis[0*vnumnodes+i];517 B_prime[(3*vnumnodes+pnumnodes)*0+3*i+1] = 0.;518 B_prime[(3*vnumnodes+pnumnodes)*0+3*i+2] = 0.;519 B_prime[(3*vnumnodes+pnumnodes)*1+3*i+0] = 0.;520 B_prime[(3*vnumnodes+pnumnodes)*1+3*i+1] = vdbasis[1*vnumnodes+i];521 B_prime[(3*vnumnodes+pnumnodes)*1+3*i+2] = 0.;522 B_prime[(3*vnumnodes+pnumnodes)*2+3*i+0] = 0.;523 B_prime[(3*vnumnodes+pnumnodes)*2+3*i+1] = 0.;524 B_prime[(3*vnumnodes+pnumnodes)*2+3*i+2] = vdbasis[2*vnumnodes+i];525 B_prime[(3*vnumnodes+pnumnodes)*3+3*i+0] = vdbasis[1*vnumnodes+i];526 B_prime[(3*vnumnodes+pnumnodes)*3+3*i+1] = vdbasis[0*vnumnodes+i];527 B_prime[(3*vnumnodes+pnumnodes)*3+3*i+2] = 0.;528 B_prime[(3*vnumnodes+pnumnodes)*4+3*i+0] = vdbasis[2*vnumnodes+i];529 B_prime[(3*vnumnodes+pnumnodes)*4+3*i+1] = 0.;530 B_prime[(3*vnumnodes+pnumnodes)*4+3*i+2] = vdbasis[0*vnumnodes+i];531 B_prime[(3*vnumnodes+pnumnodes)*5+3*i+0] = 0.;532 B_prime[(3*vnumnodes+pnumnodes)*5+3*i+1] = vdbasis[2*vnumnodes+i];533 B_prime[(3*vnumnodes+pnumnodes)*5+3*i+2] = vdbasis[1*vnumnodes+i];534 B_prime[(3*vnumnodes+pnumnodes)*6+3*i+0] = vdbasis[0*vnumnodes+i];535 B_prime[(3*vnumnodes+pnumnodes)*6+3*i+1] = vdbasis[1*vnumnodes+i];536 B_prime[(3*vnumnodes+pnumnodes)*6+3*i+2] = vdbasis[2*vnumnodes+i];537 B_prime[(3*vnumnodes+pnumnodes)*7+3*i+0] = 0.;538 B_prime[(3*vnumnodes+pnumnodes)*7+3*i+1] = 0.;539 B_prime[(3*vnumnodes+pnumnodes)*7+3*i+2] = 0.;540 }541 for(int i=0;i<pnumnodes;i++){542 B_prime[(3*vnumnodes+pnumnodes)*0+(3*vnumnodes)+i] = 0.;543 B_prime[(3*vnumnodes+pnumnodes)*1+(3*vnumnodes)+i] = 0.;544 B_prime[(3*vnumnodes+pnumnodes)*2+(3*vnumnodes)+i] = 0.;545 B_prime[(3*vnumnodes+pnumnodes)*3+(3*vnumnodes)+i] = 0.;546 B_prime[(3*vnumnodes+pnumnodes)*4+(3*vnumnodes)+i] = 0.;547 B_prime[(3*vnumnodes+pnumnodes)*5+(3*vnumnodes)+i] = 0.;548 B_prime[(3*vnumnodes+pnumnodes)*6+(3*vnumnodes)+i] = 0.;549 B_prime[(3*vnumnodes+pnumnodes)*7+(3*vnumnodes)+i] = pbasis[i];550 }551 552 /*Clean up*/553 xDelete<IssmDouble>(vdbasis);554 xDelete<IssmDouble>(pbasis);555 }556 /*}}}*/557 /*FUNCTION PentaRef::GetBprimeFSGLS {{{*/558 void PentaRef::GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss){559 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.560 * For node i, Bi' can be expressed in the actual coordinate system561 * by:562 * Bi'=[ dh/dx 0 0 0]563 * [ 0 dh/dy 0 0]564 * [ 0 0 dh/dz 0]565 * [ dh/dy dh/dx 0 0]566 * [ dh/dz 0 dh/dx 0]567 * [ 0 dh/dz dh/dy 0]568 * [ dh/dx dh/dy dh/dz 0]569 * [ 0 0 0 h]570 * where h is the interpolation function for node i.571 *572 * Same thing for the bubble fonction except that there is no fourth column573 */574 575 int i;576 IssmDouble dbasis[3][NUMNODESP1];577 IssmDouble basis[NUMNODESP1];578 579 /*Get dbasismini in actual coordinate system: */580 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);581 GetNodalFunctionsP1(basis, gauss);582 583 /*B_primeuild B_prime: */584 for(i=0;i<NUMNODESP1;i++){585 B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];586 B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;587 B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;588 B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;589 B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];590 B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;591 B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;592 B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;593 B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];594 B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = dbasis[1][i];595 B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = dbasis[0][i];596 B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;597 B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = dbasis[2][i];598 B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;599 B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = dbasis[0][i];600 B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;601 B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = dbasis[2][i];602 B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = dbasis[1][i];603 B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = dbasis[0][i];604 B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = dbasis[1][i];605 B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = dbasis[2][i];606 B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = 0.;607 B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = 0.;608 B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = 0.;609 }610 611 for(i=0;i<NUMNODESP1;i++){ //last column612 B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;613 B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;614 B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;615 B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;616 B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;617 B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;618 B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = 0.;619 B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = - basis[i];620 }621 622 }623 /*}}}*/624 /*FUNCTION PentaRef::GetBAdvec{{{*/625 void PentaRef::GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, Gauss* gauss){626 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.627 * For node i, Bi' can be expressed in the actual coordinate system628 * by:629 * Bi_advec =[ h ]630 * [ h ]631 * [ h ]632 * where h is the interpolation function for node i.633 *634 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)635 */636 637 /*Fetch number of nodes for this finite element*/638 int numnodes = this->NumberofNodes();639 640 /*Get nodal functions derivatives*/641 IssmDouble* basis=xNew<IssmDouble>(numnodes);642 GetNodalFunctions(basis,gauss);643 644 /*Build B: */645 for(int i=0;i<numnodes;i++){646 B_advec[numnodes*0+i] = basis[i];647 B_advec[numnodes*1+i] = basis[i];648 B_advec[numnodes*2+i] = basis[i];649 }650 651 /*Clean-up*/652 xDelete<IssmDouble>(basis);653 }654 /*}}}*/655 /*FUNCTION PentaRef::GetBConduct{{{*/656 void PentaRef::GetBConduct(IssmDouble* B_conduct, IssmDouble* xyz_list, Gauss* gauss){657 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.658 * For node i, Bi' can be expressed in the actual coordinate system659 * by:660 * Bi_conduct=[ dh/dx ]661 * [ dh/dy ]662 * [ dh/dz ]663 * where h is the interpolation function for node i.664 *665 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)666 */667 668 /*Fetch number of nodes for this finite element*/669 int numnodes = this->NumberofNodes();670 671 /*Get nodal functions derivatives*/672 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);673 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);674 675 /*Build B: */676 for(int i=0;i<numnodes;i++){677 B_conduct[numnodes*0+i] = dbasis[0*numnodes+i];678 B_conduct[numnodes*1+i] = dbasis[1*numnodes+i];679 B_conduct[numnodes*2+i] = dbasis[2*numnodes+i];680 }681 682 /*Clean-up*/683 xDelete<IssmDouble>(dbasis);684 }685 /*}}}*/686 /*FUNCTION PentaRef::GetBVert{{{*/687 void PentaRef::GetBVert(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){688 /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];689 where hi is the interpolation function for node i.*/690 691 /*Fetch number of nodes for this finite element*/692 int numnodes = this->NumberofNodes();693 694 /*Get nodal functions derivatives*/695 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);696 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);697 698 /*Build B: */699 for(int i=0;i<numnodes;i++){700 B[i] = dbasis[2*numnodes+i];701 }702 703 /*Clean-up*/704 xDelete<IssmDouble>(dbasis);705 }706 /*}}}*/707 /*FUNCTION PentaRef::GetBprimeAdvec{{{*/708 void PentaRef::GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, Gauss* gauss){709 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.710 * For node i, Bi' can be expressed in the actual coordinate system711 * by:712 * Biprime_advec=[ dh/dx ]713 * [ dh/dy ]714 * [ dh/dz ]715 * where h is the interpolation function for node i.716 *717 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)718 */719 720 /*Fetch number of nodes for this finite element*/721 int numnodes = this->NumberofNodes();722 723 /*Get nodal functions derivatives*/724 IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);725 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);726 727 /*Build B': */728 for(int i=0;i<numnodes;i++){729 Bprime_advec[numnodes*0+i] = dbasis[0*numnodes+i];730 Bprime_advec[numnodes*1+i] = dbasis[1*numnodes+i];731 Bprime_advec[numnodes*2+i] = dbasis[2*numnodes+i];732 }733 734 /*Clean-up*/735 xDelete<IssmDouble>(dbasis);736 }737 /*}}}*/738 /*FUNCTION PentaRef::GetBprimeVert{{{*/739 void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){740 741 GetNodalFunctions(B,gauss);742 743 }744 /*}}}*/745 /*FUNCTION PentaRef::GetBHOFriction{{{*/746 void PentaRef::GetBHOFriction(IssmDouble* B, Gauss* gauss){747 /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2x2.748 ** For node i, Bi can be expressed in the actual coordinate system749 ** by:750 ** Bi=[ N 0 ]751 ** [ 0 N ]752 ** where N is the interpolation function for node i.753 **754 ** We assume B has been allocated already, of size: 2 (2 x numnodes)755 **/756 757 /*Fetch number of nodes for this finite element*/758 int numnodes = this->NumberofNodes();759 760 /*Get nodal functions derivatives*/761 IssmDouble* basis=xNew<IssmDouble>(numnodes);762 GetNodalFunctions(basis,gauss);763 764 for(int i=0;i<numnodes;i++){765 B[2*numnodes*0+2*i+0] = basis[i];766 B[2*numnodes*0+2*i+1] = 0.;767 B[2*numnodes*1+2*i+0] = 0.;768 B[2*numnodes*1+2*i+1] = basis[i];769 }770 771 /*Clean-up*/772 xDelete<IssmDouble>(basis);773 }774 /*}}}*/775 /*FUNCTION PentaRef::GetLFS{{{*/776 void PentaRef::GetLFS(IssmDouble* LFS, Gauss* gauss){777 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.778 * For node i, Li can be expressed in the actual coordinate system779 * by:780 * Li=[ h 0 0 0 ]781 * [ 0 h 0 0 ]782 * where h is the interpolation function for node i.783 */784 785 /*Fetch number of nodes for this finite element*/786 int pnumnodes = this->NumberofNodesPressure();787 int vnumnodes = this->NumberofNodesVelocity();788 int pnumdof = pnumnodes;789 int vnumdof = vnumnodes*NDOF3;790 791 /*Get nodal functions derivatives*/792 IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);793 GetNodalFunctionsVelocity(vbasis,gauss);794 795 /*Build LFS: */796 for(int i=0;i<vnumnodes;i++){797 LFS[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];798 LFS[(vnumdof+pnumdof)*0+3*i+1] = 0.;799 LFS[(vnumdof+pnumdof)*0+3*i+2] = 0.;800 801 LFS[(vnumdof+pnumdof)*1+3*i+0] = 0.;802 LFS[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];803 LFS[(vnumdof+pnumdof)*1+3*i+2] = 0.;804 }805 806 for(int i=0;i<pnumnodes;i++){807 LFS[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;808 LFS[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;809 }810 811 /*Clean-up*/812 xDelete<IssmDouble>(vbasis);813 }814 /*}}}*/815 /*FUNCTION PentaRef::GetLprimeFS {{{*/816 void PentaRef::GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, Gauss* gauss_in){817 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.818 * For node i, Lpi can be expressed in the actual coordinate system819 * by:820 * Lpi=[ h 0 0 0]1821 * [ 0 h 0 0]2822 * [ h 0 0 0]3823 * [ 0 h 0 0]4824 * [ 0 0 h 0]5825 * [ 0 0 h 0]6826 * [ 0 0 dh/dz 0]7827 * [ 0 0 dh/dz 0]8828 * [ 0 0 dh/dz 0]9829 * [dh/dz 0 dh/dx 0]0830 * [ 0 dh/dz dh/dy 0]1831 * [ 0 0 0 h]2832 * [ 0 0 0 h]3833 * [ 0 0 0 h]4834 *835 * Li=[ h 0 0 0]1836 * [ 0 h 0 0]2837 * [ 0 0 h 0]3838 * [ 0 0 h 0]4839 * [ h 0 0 0]5840 * [ 0 h 0 0]6841 * [ h 0 0 0]7842 * [ 0 h 0 0]8843 * [ 0 0 h 0]9844 * [ 0 0 h 0]0845 * [ 0 0 h 0]1846 * [ h 0 0 0]2847 * [ 0 h 0 0]3848 * [ 0 0 h 0]4849 * where h is the interpolation function for node i.850 */851 852 int num_dof=4;853 IssmDouble L1L2l3[NUMNODESP1_2d];854 IssmDouble dbasis[3][NUMNODESP1];855 856 /*Cast gauss to GaussPenta*/857 _assert_(gauss_in->Enum()==GaussPentaEnum);858 GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);859 860 /*Get L1L2l3 in actual coordinate system: */861 L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;862 L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;863 L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;864 865 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);866 867 /*Build LprimeFS: */868 for(int i=0;i<3;i++){869 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];870 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;871 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;872 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;873 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;874 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];875 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;876 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;877 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];878 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;879 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = 0.;880 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;881 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;882 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];883 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = 0.;884 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;885 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;886 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;887 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = L1L2l3[i];888 LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;889 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;890 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;891 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = L1L2l3[i];892 LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;893 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;894 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;895 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = dbasis[2][i];896 LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = 0.;897 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;898 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;899 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = dbasis[2][i];900 LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = 0.;901 LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+0] = 0.;902 LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+1] = 0.;903 LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+2] = dbasis[2][i];904 LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+3] = 0.;905 LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+0] = dbasis[2][i];906 LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+1] = 0.;907 LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+2] = dbasis[0][i];908 LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+3] = 0.;909 LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;910 LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];911 LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];912 LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;913 LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;914 LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;915 LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;916 LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = L1L2l3[i];917 LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;918 LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;919 LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;920 LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = L1L2l3[i];921 LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;922 LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;923 LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;924 LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = L1L2l3[i];925 }926 }927 /*}}}*/928 /*FUNCTION PentaRef::GetLSSAFS {{{*/929 void PentaRef::GetLSSAFS(IssmDouble* LFS, Gauss* gauss_in){930 /*931 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.932 * For node i, Li can be expressed in the actual coordinate system933 * by:934 * Li=[ h 0 ]935 * [ 0 h ]936 * [ h 0 ]937 * [ 0 h ]938 * [ h 0 ]939 * [ 0 h ]940 * [ h 0 ]941 * [ 0 h ]942 * where h is the interpolation function for node i.943 */944 945 int num_dof=2;946 IssmDouble L1L2l3[NUMNODESP1_2d];947 948 /*Cast gauss to GaussPenta*/949 _assert_(gauss_in->Enum()==GaussPentaEnum);950 GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);951 952 /*Get L1L2l3 in actual coordinate system: */953 L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;954 L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;955 L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;956 957 /*Build LFS: */958 for(int i=0;i<3;i++){959 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];960 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;961 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;962 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];963 LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];964 LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;965 LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;966 LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];967 LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = L1L2l3[i];968 LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;969 LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;970 LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = L1L2l3[i];971 LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = L1L2l3[i];972 LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;973 LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;974 LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = L1L2l3[i];975 }976 }977 /*}}}*/978 /*FUNCTION PentaRef::GetLprimeSSAFS {{{*/979 void PentaRef::GetLprimeSSAFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, Gauss* gauss_in){980 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.981 * For node i, Lpi can be expressed in the actual coordinate system982 * by:983 * Lpi=[ h 0 0 0]984 * [ 0 h 0 0]985 * [ 0 0 h 0]986 * [ 0 0 h 0]987 * [ 0 0 dh/dz 0]988 * [ 0 0 dh/dz 0]989 * [ 0 0 0 h]990 * [ 0 0 0 h]991 * where h is the interpolation function for node i.992 */993 int num_dof=3;994 int num_dof_vel=3*NUMNODESP1b;995 int num_dof_total=3*NUMNODESP1b+1*NUMNODESP1;996 IssmDouble L1L2l3[NUMNODESP1_2d];997 IssmDouble dbasis[3][NUMNODESP1];998 999 /*Cast gauss to GaussPenta*/1000 _assert_(gauss_in->Enum()==GaussPentaEnum);1001 GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);1002 1003 /*Get L1L2l3 in actual coordinate system: */1004 L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;1005 L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;1006 L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;1007 1008 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);1009 1010 /*Build LprimeFS: */1011 for(int i=0;i<3;i++){1012 LprimeFS[num_dof_total*0+num_dof*i+0] = L1L2l3[i];1013 LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;1014 LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;1015 LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;1016 LprimeFS[num_dof_total*1+num_dof*i+1] = L1L2l3[i];1017 LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;1018 LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;1019 LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;1020 LprimeFS[num_dof_total*2+num_dof*i+2] = L1L2l3[i];1021 LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;1022 LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;1023 LprimeFS[num_dof_total*3+num_dof*i+2] = L1L2l3[i];1024 LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;1025 LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;1026 LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i];1027 LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;1028 LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;1029 LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i];1030 LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;1031 LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;1032 LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;1033 LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;1034 LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;1035 LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;1036 }1037 for(int i=3;i<7;i++){1038 LprimeFS[num_dof_total*0+num_dof*i+0] = 0.;1039 LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;1040 LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;1041 LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;1042 LprimeFS[num_dof_total*1+num_dof*i+1] = 0.;1043 LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;1044 LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;1045 LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;1046 LprimeFS[num_dof_total*2+num_dof*i+2] = 0.;1047 LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;1048 LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;1049 LprimeFS[num_dof_total*3+num_dof*i+2] = 0.;1050 LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;1051 LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;1052 LprimeFS[num_dof_total*4+num_dof*i+2] = 0.;1053 LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;1054 LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;1055 LprimeFS[num_dof_total*5+num_dof*i+2] = 0.;1056 LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;1057 LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;1058 LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;1059 LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;1060 LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;1061 LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;1062 }1063 for(int i=0;i<3;i++){1064 LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;1065 LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;1066 LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;1067 LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;1068 LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;1069 LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;1070 LprimeFS[num_dof_total*6+num_dof_vel+i] = L1L2l3[i];1071 LprimeFS[num_dof_total*7+num_dof_vel+i] = L1L2l3[i];1072 }1073 for(int i=3;i<6;i++){1074 LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;1075 LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;1076 LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;1077 LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;1078 LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;1079 LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;1080 LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.;1081 LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.;1082 }1083 }1084 /*}}}*/1085 /*FUNCTION PentaRef::GetLFSSSA {{{*/1086 void PentaRef::GetLFSSSA(IssmDouble* LFS, Gauss* gauss_in){1087 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.1088 * For node i, Li can be expressed in the actual coordinate system1089 * by:1090 * Li=[ h 0 0 ]1091 * [ 0 h 0 ]1092 * [ 0 0 h ]1093 * [ 0 0 h ]1094 * where h is the interpolation function for node i.1095 */1096 1097 int num_dof=3;1098 IssmDouble L1L2l3[NUMNODESP1_2d];1099 1100 /*Cast gauss to GaussPenta*/1101 _assert_(gauss_in->Enum()==GaussPentaEnum);1102 GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);1103 1104 /*Get L1L2l3 in actual coordinate system: */1105 L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;1106 L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;1107 L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;1108 1109 /*Build LFS: */1110 for(int i=0;i<3;i++){1111 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];1112 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;1113 LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;1114 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;1115 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];1116 LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;1117 LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;1118 LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;1119 LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];1120 LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;1121 LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;1122 LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];1123 }1124 }1125 /*}}}*/1126 /*FUNCTION PentaRef::GetLprimeFSSSA {{{*/1127 void PentaRef::GetLprimeFSSSA(IssmDouble* LprimeFS, IssmDouble* xyz_list, Gauss* gauss_in){1128 /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.1129 * For node i, Lpi can be expressed in the actual coordinate system1130 * by:1131 * Lpi=[ h 0 ]1132 * [ 0 h ]1133 * [ h 0 ]1134 * [ 0 h ]1135 * where h is the interpolation function for node i.1136 */1137 int num_dof=2;1138 IssmDouble L1L2l3[NUMNODESP1_2d];1139 IssmDouble dbasis[3][NUMNODESP1];1140 1141 /*Cast gauss to GaussPenta*/1142 _assert_(gauss_in->Enum()==GaussPentaEnum);1143 GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);1144 1145 /*Get L1L2l3 in actual coordinate system: */1146 L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;1147 L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;1148 L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;1149 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);1150 1151 /*Build LprimeFS: */1152 for(int i=0;i<3;i++){1153 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];1154 LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;1155 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;1156 LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];1157 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];1158 LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;1159 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;1160 LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];1161 }1162 }1163 /*}}}*/1164 58 /*FUNCTION PentaRef::GetJacobian {{{*/ 1165 59 void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss_in){ -
issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
r16838 r17095 41 41 void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); 42 42 void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,Gauss* gauss); 43 void GetBSSAHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);44 void GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);45 void GetBHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);46 void GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);47 void GetBFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);48 void GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);49 void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);50 void GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);51 void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);52 void GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);53 void GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);54 void GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, Gauss* gauss);55 void GetBConduct(IssmDouble* B_conduct, IssmDouble* xyz_list, Gauss* gauss);56 void GetBVert(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);57 void GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, Gauss* gauss);58 void GetBHOFriction(IssmDouble* L, Gauss* gauss);59 void GetLFS(IssmDouble* LFS, Gauss* gauss);60 void GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, Gauss* gauss);61 void GetLSSAFS(IssmDouble* LSSAFS, Gauss* gauss);62 void GetLprimeSSAFS(IssmDouble* LprimeSSAFS, IssmDouble* xyz_list, Gauss* gauss);63 void GetLFSSSA(IssmDouble* LFSSSA, Gauss* gauss);64 43 void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss); 65 44 void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss); -
issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
r16818 r17095 51 51 52 52 /*Reference Element numerics*/ 53 /*FUNCTION TriaRef::GetBHydro {{{*/54 void TriaRef::GetBHydro(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){55 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.56 * For node i, Bi can be expressed in the actual coordinate system57 * by:58 * Bi=[ dN/dx ]59 * [ dN/dy ]60 * where N is the finiteelement function for node i.61 *62 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)63 */64 65 /*Fetch number of nodes for this finite element*/66 int numnodes = this->NumberofNodes();67 68 /*Get nodal functions derivatives*/69 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);70 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);71 72 /*Build B: */73 for(int i=0;i<numnodes;i++){74 B[numnodes*0+i] = dbasis[0*numnodes+i];75 B[numnodes*1+i] = dbasis[1*numnodes+i];76 }77 78 /*Clean-up*/79 xDelete<IssmDouble>(dbasis);80 }81 /*}}}*/82 /*FUNCTION TriaRef::GetBSSA {{{*/83 void TriaRef::GetBSSA(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){84 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.85 * For node i, Bi can be expressed in the actual coordinate system86 * by:87 * Bi=[ dN/dx 0 ]88 * [ 0 dN/dy ]89 * [ 1/2*dN/dy 1/2*dN/dx ]90 * where N is the finiteelement function for node i.91 *92 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)93 */94 95 /*Fetch number of nodes for this finite element*/96 int numnodes = this->NumberofNodes();97 98 /*Get nodal functions derivatives*/99 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);100 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);101 102 /*Build B: */103 for(int i=0;i<numnodes;i++){104 B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];105 B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;106 B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;107 B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];108 B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];109 B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];110 }111 112 /*Clean-up*/113 xDelete<IssmDouble>(dbasis);114 }115 /*}}}*/116 /*FUNCTION TriaRef::GetBSSAFS {{{*/117 void TriaRef::GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){118 119 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.120 * For node i, Bi can be expressed in the actual coordinate system121 * by:122 * Bi=[ dN/dx 0 ]123 * [ 0 dN/dy ]124 * [ 1/2*dN/dy 1/2*dN/dx ]125 * where N is the finiteelement function for node i.126 *127 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)128 */129 130 /*Fetch number of nodes for this finite element*/131 int numnodes = this->NumberofNodes();132 133 /*Get nodal functions derivatives*/134 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);135 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);136 137 /*Build B': */138 for(int i=0;i<numnodes;i++){139 B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];140 B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;141 B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;142 B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];143 B[NDOF2*numnodes*2+NDOF2*i+0] = 0.5*dbasis[1*numnodes+i];144 B[NDOF2*numnodes*2+NDOF2*i+1] = 0.5*dbasis[0*numnodes+i];145 }146 147 /*Clean-up*/148 xDelete<IssmDouble>(dbasis);149 }150 /*}}}*/151 53 /*FUNCTION TriaRef::GetSegmentBFlux{{{*/ 152 54 void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2){ … … 199 101 /*Clean-up*/ 200 102 xDelete<IssmDouble>(basis); 201 }202 /*}}}*/203 /*FUNCTION TriaRef::GetBExtrusion{{{*/204 void TriaRef::GetBExtrusion(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){205 /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];206 where hi is the interpolation function for node i.*/207 208 /*Fetch number of nodes for this finite element*/209 int numnodes = this->NumberofNodes();210 211 /*Get nodal functions derivatives*/212 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);213 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);214 215 /*Build B: */216 for(int i=0;i<numnodes;i++){217 B[i] = dbasis[1*numnodes+i];218 }219 220 /*Clean-up*/221 xDelete<IssmDouble>(dbasis);222 }223 /*}}}*/224 /*FUNCTION TriaRef::GetBFS {{{*/225 void TriaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){226 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.227 * For node i, Bvi can be expressed in the actual coordinate system228 * by: Bvi=[ dphi/dx 0 ]229 * [ 0 dphi/dy ]230 * [ 1/2*dphi/dy 1/2*dphi/dx]231 * [ 0 0 ]232 * [ dphi/dx dphi/dy ]233 *234 * by: Bpi=[ 0 ]235 * [ 0 ]236 * [ 0 ]237 * [ phi_p ]238 * [ 0 ]239 * where phi is the finiteelement function for node i.240 * Same thing for Bb except the last column that does not exist.241 */242 243 /*Fetch number of nodes for this finite element*/244 int pnumnodes = this->NumberofNodesPressure();245 int vnumnodes = this->NumberofNodesVelocity();246 247 /*Get nodal functions derivatives*/248 IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);249 IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);250 GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);251 GetNodalFunctionsPressure(pbasis,gauss);252 253 /*Build B: */254 for(int i=0;i<vnumnodes;i++){255 B[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];256 B[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;257 B[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;258 B[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];259 B[(2*vnumnodes+pnumnodes)*2+2*i+0] = .5*vdbasis[1*vnumnodes+i];260 B[(2*vnumnodes+pnumnodes)*2+2*i+1] = .5*vdbasis[0*vnumnodes+i];261 B[(2*vnumnodes+pnumnodes)*3+2*i+0] = 0.;262 B[(2*vnumnodes+pnumnodes)*3+2*i+1] = 0.;263 B[(2*vnumnodes+pnumnodes)*4+2*i+0] = vdbasis[0*vnumnodes+i];264 B[(2*vnumnodes+pnumnodes)*4+2*i+1] = vdbasis[1*vnumnodes+i];265 }266 for(int i=0;i<pnumnodes;i++){267 B[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;268 B[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;269 B[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;270 B[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = pbasis[i];271 B[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = 0.;272 }273 274 /*Clean up*/275 xDelete<IssmDouble>(vdbasis);276 xDelete<IssmDouble>(pbasis);277 }278 /*}}}*/279 /*FUNCTION TriaRef::GetBprimeFS {{{*/280 void TriaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss){281 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.282 * For node i, Bi' can be expressed in the actual coordinate system283 * by:284 * Bvi' = [ dphi/dx 0 ]285 * [ 0 dphi/dy ]286 * [ dphi/dy dphi/dx ]287 * [ dphi/dx dphi/dy ]288 * [ 0 0 ]289 *290 * by: Bpi=[ 0 ]291 * [ 0 ]292 * [ 0 ]293 * [ 0 ]294 * [ phi ]295 * where phi is the finiteelement function for node i.296 */297 298 /*Fetch number of nodes for this finite element*/299 int pnumnodes = this->NumberofNodesPressure();300 int vnumnodes = this->NumberofNodesVelocity();301 302 /*Get nodal functions derivatives*/303 IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);304 IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);305 GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);306 GetNodalFunctionsPressure(pbasis,gauss);307 308 /*Build B_prime: */309 for(int i=0;i<vnumnodes;i++){310 B_prime[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];311 B_prime[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;312 B_prime[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;313 B_prime[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];314 B_prime[(2*vnumnodes+pnumnodes)*2+2*i+0] = vdbasis[1*vnumnodes+i];315 B_prime[(2*vnumnodes+pnumnodes)*2+2*i+1] = vdbasis[0*vnumnodes+i];316 B_prime[(2*vnumnodes+pnumnodes)*3+2*i+0] = vdbasis[0*vnumnodes+i];317 B_prime[(2*vnumnodes+pnumnodes)*3+2*i+1] = vdbasis[1*vnumnodes+i];318 B_prime[(2*vnumnodes+pnumnodes)*4+2*i+0] = 0.;319 B_prime[(2*vnumnodes+pnumnodes)*4+2*i+1] = 0.;320 }321 for(int i=0;i<pnumnodes;i++){322 B_prime[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;323 B_prime[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;324 B_prime[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;325 B_prime[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = 0.;326 B_prime[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = pbasis[i];327 }328 329 /*Clean up*/330 xDelete<IssmDouble>(vdbasis);331 xDelete<IssmDouble>(pbasis);332 }333 /*}}}*/334 /*FUNCTION TriaRef::GetBMasstransport{{{*/335 void TriaRef::GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){336 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.337 * For node i, Bi can be expressed in the actual coordinate system338 * by:339 * Bi=[ N ]340 * [ N ]341 * where N is the finiteelement function for node i.342 *343 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)344 */345 346 /*Fetch number of nodes for this finite element*/347 int numnodes = this->NumberofNodes();348 349 /*Get nodal functions*/350 IssmDouble* basis=xNew<IssmDouble>(numnodes);351 GetNodalFunctions(basis,gauss);352 353 /*Build B: */354 for(int i=0;i<numnodes;i++){355 B[numnodes*0+i] = basis[i];356 B[numnodes*1+i] = basis[i];357 }358 359 /*Clean-up*/360 xDelete<IssmDouble>(basis);361 }362 /*}}}*/363 /*FUNCTION TriaRef::GetBprimeSSA {{{*/364 void TriaRef::GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){365 366 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.367 * For node i, Bi' can be expressed in the actual coordinate system368 * by:369 * Bi_prime=[ 2*dN/dx dN/dy ]370 * [ dN/dx 2*dN/dy ]371 * [ dN/dy dN/dx ]372 * where hNis the finiteelement function for node i.373 *374 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)375 */376 377 /*Fetch number of nodes for this finite element*/378 int numnodes = this->NumberofNodes();379 380 /*Get nodal functions derivatives*/381 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);382 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);383 384 /*Build B': */385 for(int i=0;i<numnodes;i++){386 Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];387 Bprime[NDOF2*numnodes*0+NDOF2*i+1] = dbasis[1*numnodes+i];388 Bprime[NDOF2*numnodes*1+NDOF2*i+0] = dbasis[0*numnodes+i];389 Bprime[NDOF2*numnodes*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];390 Bprime[NDOF2*numnodes*2+NDOF2*i+0] = dbasis[1*numnodes+i];391 Bprime[NDOF2*numnodes*2+NDOF2*i+1] = dbasis[0*numnodes+i];392 }393 394 /*Clean-up*/395 xDelete<IssmDouble>(dbasis);396 }397 /*}}}*/398 /*FUNCTION TriaRef::GetBprimeSSAFS {{{*/399 void TriaRef::GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){400 /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.401 * For node i, Bprimei can be expressed in the actual coordinate system402 * by:403 * Bprimei=[ dN/dx 0 ]404 * [ 0 dN/dy ]405 * [ dN/dy dN/dx ]406 N [ dN/dx dN/dy ]407 * where N is the finiteelement function for node i.408 *409 * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)410 */411 412 /*Fetch number of nodes for this finite element*/413 int numnodes = this->NumberofNodes();414 415 /*Get nodal functions*/416 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);417 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);418 419 /*Build Bprime: */420 for(int i=0;i<numnodes;i++){421 Bprime[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];422 Bprime[NDOF2*numnodes*0+NDOF2*i+1] = 0.;423 Bprime[NDOF2*numnodes*1+NDOF2*i+0] = 0.;424 Bprime[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];425 Bprime[NDOF2*numnodes*2+NDOF2*i+0] = dbasis[1*numnodes+i];426 Bprime[NDOF2*numnodes*2+NDOF2*i+1] = dbasis[0*numnodes+i];427 Bprime[NDOF2*numnodes*3+NDOF2*i+0] = dbasis[0*numnodes+i];428 Bprime[NDOF2*numnodes*3+NDOF2*i+1] = dbasis[1*numnodes+i];429 }430 431 /*Clean-up*/432 xDelete<IssmDouble>(dbasis);433 }434 /*}}}*/435 /*FUNCTION TriaRef::GetBprimeMasstransport{{{*/436 void TriaRef::GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){437 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.438 * For node i, Bi' can be expressed in the actual coordinate system439 * by:440 * Bi_prime=[ dN/dx ]441 * [ dN/dy ]442 * where N is the finiteelement function for node i.443 *444 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)445 */446 447 /*Fetch number of nodes for this finite element*/448 int numnodes = this->NumberofNodes();449 450 /*Get nodal functions derivatives*/451 IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);452 GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);453 454 /*Build B': */455 for(int i=0;i<numnodes;i++){456 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];457 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];458 }459 460 /*Clean-up*/461 xDelete<IssmDouble>(dbasis);462 }463 /*}}}*/464 /*FUNCTION TriaRef::GetBSSAFriction{{{*/465 void TriaRef::GetBSSAFriction(IssmDouble* B, IssmDouble* xyz_list,Gauss* gauss){466 /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2.467 * For node i, Bi can be expressed in the actual coordinate system468 * by:469 * Bi=[ N 0 ]470 * [ 0 N ]471 * where N is the finiteelement function for node i.472 *473 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)474 */475 476 /*Fetch number of nodes for this finite element*/477 int numnodes = this->NumberofNodes();478 479 /*Get nodal functions derivatives*/480 IssmDouble* basis=xNew<IssmDouble>(numnodes);481 GetNodalFunctions(basis,gauss);482 483 /*Build L: */484 for(int i=0;i<numnodes;i++){485 B[2*numnodes*0+2*i+0] = basis[i];486 B[2*numnodes*0+2*i+1] = 0.;487 B[2*numnodes*1+2*i+0] = 0.;488 B[2*numnodes*1+2*i+1] = basis[i];489 }490 491 /*Clean-up*/492 xDelete<IssmDouble>(basis);493 }494 /*}}}*/495 /*FUNCTION TriaRef::GetLFS{{{*/496 void TriaRef::GetLFS(IssmDouble* LFS, Gauss* gauss){497 /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.498 * For node i, Li can be expressed in the actual coordinate system499 * by:500 * Li=[ h 0 0 ]501 * where h is the interpolation function for node i.502 */503 504 /*Fetch number of nodes for this finite element*/505 int pnumnodes = this->NumberofNodesPressure();506 int vnumnodes = this->NumberofNodesVelocity();507 int pnumdof = pnumnodes;508 int vnumdof = vnumnodes*NDOF2;509 510 /*Get nodal functions derivatives*/511 IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);512 GetNodalFunctionsVelocity(vbasis,gauss);513 514 /*Build LFS: */515 for(int i=0;i<vnumnodes;i++){516 LFS[2*i+0] = vbasis[i];517 LFS[2*i+1] = 0.;518 }519 520 for(int i=0;i<pnumnodes;i++){521 LFS[i+vnumdof+0] = 0.;522 }523 524 /*Clean-up*/525 xDelete<IssmDouble>(vbasis);526 103 } 527 104 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
r16818 r17095 23 23 24 24 /*Numerics*/ 25 void GetBExtrusion(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);26 void GetBFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);27 void GetBSSA(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);28 void GetBSSAFS(IssmDouble* B , IssmDouble* xyz_list, Gauss* gauss);29 void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);30 void GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);31 void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);32 void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);33 void GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);34 void GetBHydro(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);35 void GetBSSAFriction(IssmDouble* L, IssmDouble* xyz_list,Gauss* gauss);36 25 void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss); 37 26 void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); 38 27 void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); 39 28 void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,Gauss* gauss); 40 void GetLFS(IssmDouble* LFS, Gauss* gauss);41 29 void GetNodalFunctions(IssmDouble* basis,Gauss* gauss); 42 30 void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);
Note:
See TracChangeset
for help on using the changeset viewer.
