Changeset 5211
- Timestamp:
- 08/12/10 15:57:40 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r5207 r5211 3321 3321 D_scalar=D_scalar*dt; 3322 3322 } 3323 K[0][0]=D_scalar*pow(u,2); K[0][1]=D_scalar*fabs(u)*fabs(v); 3324 K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2); 3323 K[0][0]=D_scalar*pow(u,2); K[0][1]=D_scalar*(u)*(v); 3324 K[1][0]=D_scalar*(u)*(v); K[1][1]=D_scalar*pow(v,2); 3325 //K[0][0]=D_scalar*pow(u,2); K[0][1]=D_scalar*fabs(u)*fabs(v); 3326 //K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2); 3325 3327 3326 3328 /*Get B_artdiff: */ … … 5191 5193 } 5192 5194 } 5193 else {5195 else if (approximation==PattynApproximationEnum){ 5194 5196 InputUpdateFromSolutionDiagnosticPattyn(solution); 5197 } 5198 else if (approximation==MacAyealPattynApproximationEnum){ 5199 InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution); 5195 5200 } 5196 5201 } … … 5288 5293 /*Free ressources:*/ 5289 5294 xfree((void**)&doflist); 5295 } 5296 /*}}}*/ 5297 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{1*/ 5298 void Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(double* solution){ 5299 5300 int i; 5301 5302 const int numvertices=6; 5303 const int numdofpervertex=4; 5304 const int solutionnumdof=2; 5305 const int numdof=solutionnumdof*numvertices; 5306 int* doflist=NULL; 5307 int* doflistbase=NULL; 5308 double macayeal_values[numdof]; 5309 double pattyn_values[numdof]; 5310 double vx[numvertices]; 5311 double vy[numvertices]; 5312 double vz[numvertices]; 5313 double vel[numvertices]; 5314 int dummy; 5315 double pressure[numvertices]; 5316 double surface[numvertices]; 5317 double rho_ice,g; 5318 double xyz_list[numvertices][3]; 5319 double gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 5320 5321 Input *vz_input = NULL; 5322 double *vz_ptr = NULL; 5323 Penta *penta = NULL; 5324 5325 /*OK, we have to add results of this element for pattyn 5326 * and results from the penta at base for macayeal. Now recover results*/ 5327 penta=GetBasalElement(); 5328 5329 /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */ 5330 GetDofList(&doflist); 5331 penta->GetDofList(&doflistbase); 5332 5333 /*Get node data: */ 5334 GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); 5335 5336 /*Use the dof list to index into the solution vector: */ 5337 for(i=0;i<numdof;i++){ 5338 pattyn_values[i]=solution[doflist[i]]; 5339 macayeal_values[i]=solution[doflistbase[i]]; 5340 } 5341 5342 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ 5343 for(i=0;i<numvertices;i++){ 5344 vx[i]=macayeal_values[i*numdofpervertex+0]+pattyn_values[i*numdofpervertex+0]; 5345 vy[i]=macayeal_values[i*numdofpervertex+1]+pattyn_values[i*numdofpervertex+1]; 5346 } 5347 5348 /*Get Vz*/ 5349 vz_input=inputs->GetInput(VzEnum); 5350 if (vz_input){ 5351 if (vz_input->Enum()!=PentaVertexInputEnum){ 5352 ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum())); 5353 } 5354 vz_input->GetValuesPtr(&vz_ptr,&dummy); 5355 for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i]; 5356 } 5357 else{ 5358 for(i=0;i<numvertices;i++) vz[i]=0.0; 5359 } 5360 5361 /*Now Compute vel*/ 5362 for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5); 5363 5364 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 5365 *so the pressure is just the pressure at the z elevation: */ 5366 rho_ice=matpar->GetRhoIce(); 5367 g=matpar->GetG(); 5368 inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); 5369 5370 for(i=0;i<numvertices;i++){ 5371 pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); 5372 } 5373 /*Now, we have to move the previous Vx and Vy inputs to old 5374 * status, otherwise, we'll wipe them off: */ 5375 this->inputs->ChangeEnum(VxEnum,VxOldEnum); 5376 this->inputs->ChangeEnum(VyEnum,VyOldEnum); 5377 this->inputs->ChangeEnum(PressureEnum,PressureOldEnum); 5378 5379 /*Add vx and vy as inputs to the tria element: */ 5380 this->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); 5381 this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); 5382 this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); 5383 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 5384 5385 /*Free ressources:*/ 5386 xfree((void**)&doflist); 5387 xfree((void**)&doflistbase); 5290 5388 } 5291 5389 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Penta.h
r5203 r5211 167 167 void InputUpdateFromSolutionDiagnosticHoriz( double* solutiong); 168 168 void InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong); 169 void InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong); 169 170 void InputUpdateFromSolutionDiagnosticPattyn( double* solutiong); 170 171 void InputUpdateFromSolutionDiagnosticHutter( double* solutiong);
Note:
See TracChangeset
for help on using the changeset viewer.