Changeset 15486
- Timestamp:
- 07/11/13 16:30:33 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r15485 r15486 247 247 /*Fetch number of nodes and dof for this finite element*/ 248 248 int numnodes = this->NumberofNodes(); 249 int numdof = numnodes*NDOF 2;249 int numdof = numnodes*NDOF1; 250 250 251 251 /*Initialize Element matrix and vectors*/ 252 252 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); 253 IssmDouble* basis = xNew<IssmDouble>(num dof);253 IssmDouble* basis = xNew<IssmDouble>(numnodes); 254 254 255 255 /*Retrieve all inputs and parameters*/ … … 266 266 D=gauss->weight*Jdet; 267 267 268 TripleMultiply(basis,1,num nodes,1,268 TripleMultiply(basis,1,numdof,1, 269 269 &D,1,1,0, 270 basis,1,num nodes,0,270 basis,1,numdof,0, 271 271 &Ke->values[0],1); 272 272 } … … 6244 6244 6245 6245 switch(GetElementType()){ 6246 case P1Enum: 6246 case P1Enum: case P2Enum: 6247 6247 return CreateKMatrixPrognostic_CG(); 6248 6248 case P1DGEnum: … … 6256 6256 ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){ 6257 6257 6258 /*Constants*/6259 const int numdof=NDOF1*NUMVERTICES;6260 6261 6258 /*Intermediaries */ 6262 6259 int stabilization; 6263 6260 int dim; 6264 IssmDouble Jdettria,D L_scalar,dt,h;6261 IssmDouble Jdettria,D_scalar,dt,h; 6265 6262 IssmDouble vel,vx,vy,dvxdx,dvydy; 6266 6263 IssmDouble dvx[2],dvy[2]; 6267 6264 IssmDouble v_gauss[2]={0.0}; 6268 6265 IssmDouble xyz_list[NUMVERTICES][3]; 6269 IssmDouble basis[NUMVERTICES]; 6270 IssmDouble B[2][NUMVERTICES]; 6271 IssmDouble Bprime[2][NUMVERTICES]; 6272 IssmDouble K[2][2] ={0.0}; 6273 IssmDouble KDL[2][2] ={0.0}; 6274 IssmDouble DL[2][2] ={0.0}; 6275 IssmDouble DLprime[2][2] ={0.0}; 6276 6277 /*Initialize Element matrix*/ 6278 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); 6266 6267 /*Fetch number of nodes and dof for this finite element*/ 6268 int numnodes = this->NumberofNodes(); 6269 int numdof = numnodes*NDOF1; 6270 6271 /*Initialize Element matrix and vectors*/ 6272 ElementMatrix* Ke = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum); 6273 IssmDouble* basis = xNew<IssmDouble>(numnodes); 6274 IssmDouble* B = xNew<IssmDouble>(2*numdof); 6275 IssmDouble* Bprime = xNew<IssmDouble>(2*numdof); 6276 IssmDouble D[2][2]; 6279 6277 6280 6278 /*Retrieve all inputs and parameters*/ … … 6302 6300 6303 6301 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); 6304 GetNodalFunctions( &basis[0],gauss);6302 GetNodalFunctions(basis,gauss); 6305 6303 6306 6304 vxaverage_input->GetInputValue(&vx,gauss); … … 6309 6307 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); 6310 6308 6311 D L_scalar=gauss->weight*Jdettria;6312 6313 TripleMultiply( &basis[0],1,numdof,1,6314 &D L_scalar,1,1,0,6315 &basis[0],1,numdof,0,6309 D_scalar=gauss->weight*Jdettria; 6310 6311 TripleMultiply(basis,1,numdof,1, 6312 &D_scalar,1,1,0, 6313 basis,1,numdof,0, 6316 6314 &Ke->values[0],1); 6317 6315 6318 GetBPrognostic( &B[0][0], &xyz_list[0][0],gauss);6319 GetBprimePrognostic( &Bprime[0][0], &xyz_list[0][0],gauss);6316 GetBPrognostic(B,&xyz_list[0][0],gauss); 6317 GetBprimePrognostic(Bprime,&xyz_list[0][0],gauss); 6320 6318 6321 6319 dvxdx=dvx[0]; 6322 6320 dvydy=dvy[1]; 6323 DL_scalar=dt*gauss->weight*Jdettria; 6324 6325 DL[0][0]=DL_scalar*dvxdx; 6326 DL[1][1]=DL_scalar*dvydy; 6327 DLprime[0][0]=DL_scalar*vx; 6328 DLprime[1][1]=DL_scalar*vy; 6329 6330 TripleMultiply( &B[0][0],2,numdof,1, 6331 &DL[0][0],2,2,0, 6332 &B[0][0],2,numdof,0, 6321 D_scalar=dt*gauss->weight*Jdettria; 6322 6323 D[0][0]=D_scalar*dvxdx; 6324 D[0][1]=0.; 6325 D[1][1]=D_scalar*dvydy; 6326 D[1][0]=0.; 6327 TripleMultiply(B,2,numdof,1, 6328 &D[0][0],2,2,0, 6329 B,2,numdof,0, 6333 6330 &Ke->values[0],1); 6334 6331 6335 TripleMultiply( &B[0][0],2,numdof,1, 6336 &DLprime[0][0],2,2,0, 6337 &Bprime[0][0],2,numdof,0, 6332 D[0][0]=D_scalar*vx; 6333 D[1][1]=D_scalar*vy; 6334 TripleMultiply(B,2,numdof,1, 6335 &D[0][0],2,2,0, 6336 Bprime,2,numdof,0, 6338 6337 &Ke->values[0],1); 6339 6338 … … 6341 6340 /*Streamline upwinding*/ 6342 6341 vel=sqrt(vx*vx+vy*vy)+1.e-8; 6343 K[0][0]=h/(2*vel)*vx*vx;6344 K[1][0]=h/(2*vel)*vy*vx;6345 K[0][1]=h/(2*vel)*vx*vy;6346 K[1][1]=h/(2*vel)*vy*vy;6342 D[0][0]=h/(2*vel)*vx*vx; 6343 D[1][0]=h/(2*vel)*vy*vx; 6344 D[0][1]=h/(2*vel)*vx*vy; 6345 D[1][1]=h/(2*vel)*vy*vy; 6347 6346 } 6348 6347 else if(stabilization==1){ … … 6350 6349 vxaverage_input->GetInputAverage(&vx); 6351 6350 vyaverage_input->GetInputAverage(&vy); 6352 K[0][0]=h/2.0*fabs(vx);6353 K[0][1]=0.;6354 K[1][0]=0.;6355 K[1][1]=h/2.0*fabs(vy);6351 D[0][0]=h/2.0*fabs(vx); 6352 D[0][1]=0.; 6353 D[1][0]=0.; 6354 D[1][1]=h/2.0*fabs(vy); 6356 6355 } 6357 6356 if(stabilization==1 || stabilization==2){ 6358 KDL[0][0]=DL_scalar*K[0][0];6359 KDL[1][0]=DL_scalar*K[1][0];6360 KDL[0][1]=DL_scalar*K[0][1];6361 KDL[1][1]=DL_scalar*K[1][1];6362 TripleMultiply( &Bprime[0][0],2,numdof,1,6363 & KDL[0][0],2,2,0,6364 &Bprime[0][0],2,numdof,0,6357 D[0][0]=D_scalar*D[0][0]; 6358 D[1][0]=D_scalar*D[1][0]; 6359 D[0][1]=D_scalar*D[0][1]; 6360 D[1][1]=D_scalar*D[1][1]; 6361 TripleMultiply(Bprime,2,numdof,1, 6362 &D[0][0],2,2,0, 6363 Bprime,2,numdof,0, 6365 6364 &Ke->values[0],1); 6366 6365 } … … 6369 6368 /*Clean up and return*/ 6370 6369 delete gauss; 6370 xDelete<IssmDouble>(basis); 6371 xDelete<IssmDouble>(B); 6372 xDelete<IssmDouble>(Bprime); 6371 6373 return Ke; 6372 6374 } … … 6451 6453 6452 6454 switch(GetElementType()){ 6453 case P1Enum: 6455 case P1Enum: case P2Enum: 6454 6456 return CreatePVectorPrognostic_CG(); 6455 6457 case P1DGEnum: … … 6463 6465 ElementVector* Tria::CreatePVectorPrognostic_CG(void){ 6464 6466 6465 /*Constants*/6466 const int numdof=NDOF1*NUMVERTICES;6467 6468 6467 /*Intermediaries */ 6469 6468 IssmDouble Jdettria,dt; 6470 6469 IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g; 6471 6470 IssmDouble xyz_list[NUMVERTICES][3]; 6472 IssmDouble basis[NUMVERTICES]; 6473 GaussTria* gauss=NULL; 6474 6475 /*Initialize Element vector*/ 6476 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); 6471 6472 /*Fetch number of nodes and dof for this finite element*/ 6473 int numnodes = this->NumberofNodes(); 6474 int numdof = numnodes*NDOF1; 6475 6476 /*Initialize Element vector and other vectors*/ 6477 ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); 6478 IssmDouble* basis = xNew<IssmDouble>(numnodes); 6477 6479 6478 6480 /*Retrieve all inputs and parameters*/ 6479 6481 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 6480 6482 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); 6481 Input* surface_mass_balance_input =inputs->GetInput(SurfaceforcingsMassBalanceEnum);_assert_(surface_mass_balance_input);6482 Input* basal_melting_input =inputs->GetInput(BasalforcingsMeltingRateEnum);_assert_(basal_melting_input);6483 Input* basal_melting_correction_input =inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);6484 Input* thickness_input =inputs->GetInput(ThicknessEnum);_assert_(thickness_input);6483 Input* surface_mass_balance_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input); 6484 Input* basal_melting_input = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input); 6485 Input* basal_melting_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum); 6486 Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input); 6485 6487 6486 6488 /*Initialize basal_melting_correction_g to 0, do not forget!:*/ 6487 6489 /* Start looping on the number of gaussian points: */ 6488 gauss=new GaussTria(2);6490 GaussTria* gauss=new GaussTria(2); 6489 6491 for(int ig=gauss->begin();ig<gauss->end();ig++){ 6490 6492 … … 6492 6494 6493 6495 GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss); 6494 GetNodalFunctions( &basis[0],gauss);6496 GetNodalFunctions(basis,gauss); 6495 6497 6496 6498 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss); … … 6507 6509 /*Clean up and return*/ 6508 6510 delete gauss; 6511 xDelete<IssmDouble>(basis); 6509 6512 return pe; 6510 6513 }
Note:
See TracChangeset
for help on using the changeset viewer.