Changeset 5096 for issm/trunk/src/c/objects/Elements/Penta.cpp
- Timestamp:
- 08/09/10 14:48:52 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r5075 r5096 2048 2048 const int numdof=2*numgrids; 2049 2049 double xyz_list[numgrids][3]; 2050 int doflist[numdof]; 2051 int numberofdofspernode; 2050 int* doflist=NULL; 2052 2051 2053 2052 /* 3d gaussian points: */ … … 2172 2171 /* Get node coordinates and dof list: */ 2173 2172 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2174 GetDofList(&doflist [0],&numberofdofspernode);2173 GetDofList(&doflist); 2175 2174 2176 2175 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 2274 2273 xfree((void**)&third_gauss_area_coord2d); 2275 2274 xfree((void**)&gauss_weights2d); 2276 2275 xfree((void**)&doflist); 2277 2276 } 2278 2277 /*}}}*/ … … 2285 2284 const int NDOF2=2; 2286 2285 const int numdofs=NDOF2*numgrids; 2287 int doflist[numdofs];2286 int* doflist=NULL; 2288 2287 double Ke_gg[numdofs][numdofs]={0.0}; 2289 int numberofdofspernode;2290 2288 bool onbed; 2291 2289 bool onsurface; … … 2306 2304 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 2307 2305 2308 GetDofList(&doflist [0],&numberofdofspernode);2306 GetDofList(&doflist); 2309 2307 2310 2308 /*Spawn 3 beam elements: */ … … 2357 2355 MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES); 2358 2356 2357 /*Free ressources:*/ 2358 xfree((void**)&doflist); 2359 2359 2360 } 2360 2361 /*}}}*/ … … 2367 2368 const int DOFPERGRID=4; 2368 2369 const int numdof=numgrids*DOFPERGRID; 2369 int doflist[numdof]; 2370 int numberofdofspernode; 2370 int* doflist=NULL; 2371 2371 2372 2372 const int numgrids2d=3; … … 2469 2469 /* Get node coordinates and dof list: */ 2470 2470 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2471 GetDofList(&doflist [0],&numberofdofspernode);2471 GetDofList(&doflist); 2472 2472 2473 2473 /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 2636 2636 xfree((void**)&vert_gauss_coord); 2637 2637 xfree((void**)&vert_gauss_weights); 2638 xfree((void**)&doflist); 2638 2639 2639 2640 return; … … 2651 2652 const int numdof=NDOF1*numgrids; 2652 2653 double xyz_list[numgrids][3]; 2653 int doflist[numdof]; 2654 int numberofdofspernode; 2654 int* doflist=NULL; 2655 2655 2656 2656 /* 3d gaussian points: */ … … 2704 2704 /* Get node coordinates and dof list: */ 2705 2705 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2706 GetDofList(&doflist [0],&numberofdofspernode);2706 GetDofList(&doflist); 2707 2707 2708 2708 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 2763 2763 xfree((void**)&area_gauss_weights); 2764 2764 xfree((void**)&vert_gauss_weights); 2765 xfree((void**)&doflist); 2765 2766 } 2766 2767 /*}}}*/ … … 2870 2871 const int NDOF1=1; 2871 2872 const int numdof=NDOF1*numgrids; 2872 double xyz_list[numgrids][3]; 2873 int doflist[numdof]; 2874 int numberofdofspernode; 2873 double xyz_list[numgrids][3]; 2874 int* doflist=NULL; 2875 2875 2876 2876 /* gaussian points: */ … … 2947 2947 /* Get node coordinates and dof list: */ 2948 2948 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2949 GetDofList(&doflist [0],&numberofdofspernode);2949 GetDofList(&doflist); 2950 2950 2951 2951 // /*recovre material parameters: */ … … 3088 3088 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); 3089 3089 3090 //Ice/ocean heat exchange flux on ice shelf base 3091 if(onbed && shelf){ 3092 3093 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 3094 tria->CreateKMatrixThermal(Kgg); 3095 delete tria->matice; delete tria; 3096 } 3097 3098 /*Free ressources:*/ 3090 3099 xfree((void**)&first_gauss_area_coord); 3091 3100 xfree((void**)&second_gauss_area_coord); … … 3094 3103 xfree((void**)&vert_gauss_weights); 3095 3104 xfree((void**)&vert_gauss_coord); 3096 3097 //Ice/ocean heat exchange flux on ice shelf base 3098 if(onbed && shelf){ 3099 3100 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 3101 tria->CreateKMatrixThermal(Kgg); 3102 delete tria->matice; delete tria; 3103 } 3105 xfree((void**)&doflist); 3106 3104 3107 } 3105 3108 /*}}}*/ … … 3186 3189 const int numdof=NDOF2*numgrids; 3187 3190 double xyz_list[numgrids][3]; 3188 int doflist[numdof]; 3189 int numberofdofspernode; 3191 int* doflist=NULL; 3190 3192 3191 3193 /* parameters: */ … … 3262 3264 /* Get node coordinates and dof list: */ 3263 3265 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 3264 GetDofList(&doflist [0],&numberofdofspernode);3266 GetDofList(&doflist); 3265 3267 3266 3268 /*Get gaussian points and weights :*/ … … 3327 3329 xfree((void**)&area_gauss_weights); 3328 3330 xfree((void**)&vert_gauss_weights); 3329 3331 xfree((void**)&doflist); 3330 3332 } 3331 3333 /*}}}*/ … … 3383 3385 const int NDOF2=2; 3384 3386 const int numdofs=NDOF2*numgrids; 3385 int doflist[numdofs]; 3386 int numberofdofspernode; 3387 int* doflist=NULL; 3387 3388 double pe_g[numdofs]={0.0}; 3388 3389 double xyz_list[numgrids][3]; … … 3428 3429 3429 3430 /*recover doflist: */ 3430 GetDofList(&doflist [0],&numberofdofspernode);3431 GetDofList(&doflist); 3431 3432 3432 3433 /*recover some inputs: */ … … 3523 3524 xfree((void**)&gauss_weights); 3524 3525 xfree((void**)&segment_gauss_coord); 3526 xfree((void**)&doflist); 3525 3527 } 3526 3528 /*}}}*/ … … 3535 3537 const int numdof=numgrids*DOFPERGRID; 3536 3538 const int numgrids2d=3; 3537 int numdof2d=numgrids2d*DOFPERGRID; 3538 int doflist[numdof]; 3539 int numberofdofspernode; 3539 int numdof2d=numgrids2d*DOFPERGRID; 3540 int* doflist=NULL; 3540 3541 3541 3542 /*Material properties: */ … … 3627 3628 /* Get node coordinates and dof list: */ 3628 3629 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 3629 GetDofList(&doflist [0],&numberofdofspernode);3630 GetDofList(&doflist); 3630 3631 3631 3632 /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 3777 3778 xfree((void**)&area_gauss_weights); 3778 3779 xfree((void**)&vert_gauss_coord); 3779 xfree((void**)&vert_gauss_weights); 3780 xfree((void**)&vert_gauss_weights); 3781 xfree((void**)&doflist); 3780 3782 3781 3783 } … … 3815 3817 const int numdof=NDOF1*numgrids; 3816 3818 double xyz_list[numgrids][3]; 3817 int doflist[numdof]; 3818 int numberofdofspernode; 3819 int* doflist=NULL; 3819 3820 3820 3821 /* gaussian points: */ … … 3876 3877 /* Get node coordinates and dof list: */ 3877 3878 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 3878 GetDofList(&doflist [0],&numberofdofspernode);3879 GetDofList(&doflist); 3879 3880 3880 3881 /*Get gaussian points and weights :*/ … … 3935 3936 xfree((void**)&area_gauss_weights); 3936 3937 xfree((void**)&vert_gauss_weights); 3938 xfree((void**)&doflist); 3937 3939 } 3938 3940 /*}}}*/ … … 4015 4017 const int NDOF1=1; 4016 4018 const int numdof=numgrids*NDOF1; 4017 int doflist[numdof]; 4018 int numberofdofspernode; 4019 int* doflist=NULL; 4019 4020 4020 4021 /*Grid data: */ … … 4096 4097 /* Get node coordinates and dof list: */ 4097 4098 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 4098 GetDofList(&doflist [0],&numberofdofspernode);4099 GetDofList(&doflist); 4099 4100 4100 4101 /*recovre material parameters: */ … … 4191 4192 xfree((void**)&area_gauss_weights); 4192 4193 xfree((void**)&vert_gauss_weights); 4194 xfree((void**)&doflist); 4193 4195 4194 4196 } … … 4234 4236 /*}}}*/ 4235 4237 /*FUNCTION Penta::GetDofList {{{1*/ 4236 void Penta::GetDofList(int* doflist,int* pnumberofdofspernode){4238 void Penta::GetDofList(int** pdoflist){ 4237 4239 4238 4240 int i,j; 4239 int doflist_per_node[MAXDOFSPERNODE]; 4240 int numberofdofspernode; 4241 4241 int numberofdofs=0; 4242 int count=0; 4243 4244 /*output: */ 4245 int* doflist=NULL; 4246 4247 /*First, figure out size of doflist: */ 4242 4248 for(i=0;i<6;i++){ 4243 nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); 4244 for(j=0;j<numberofdofspernode;j++){ 4245 doflist[i*numberofdofspernode+j]=doflist_per_node[j]; 4246 } 4249 numberofdofs+=nodes[i]->GetNumberOfDofs(); 4250 } 4251 4252 /*Allocate: */ 4253 doflist=(int*)xmalloc(numberofdofs*sizeof(int)); 4254 4255 /*Populate: */ 4256 count=0; 4257 for(i=0;i<6;i++){ 4258 nodes[i]->GetDofList(doflist+count); 4259 count+=nodes[i]->GetNumberOfDofs(); 4247 4260 } 4248 4261 4249 4262 /*Assign output pointers:*/ 4250 *p numberofdofspernode=numberofdofspernode;4263 *pdoflist=doflist; 4251 4264 4252 4265 } … … 4440 4453 const int numdof=numdofpervertex*numvertices; 4441 4454 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}}; 4442 4443 int doflist[numdof]; 4455 int* doflist=NULL; 4444 4456 double values[numdof]; 4445 4457 double vx; 4446 4458 double vy; 4447 4459 4448 int dummy;4449 4450 4460 /*Get dof list: */ 4451 GetDofList(&doflist [0],&dummy);4461 GetDofList(&doflist); 4452 4462 4453 4463 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ … … 4465 4475 VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); 4466 4476 4477 /*Free ressources:*/ 4478 xfree((void**)&doflist); 4467 4479 } 4468 4480 /*}}}*/ … … 4476 4488 const int numdof=numdofpervertex*numvertices; 4477 4489 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}}; 4478 4479 int doflist[numdof]; 4490 int* doflist=NULL; 4480 4491 double values[numdof]; 4481 4492 double vx; 4482 4493 double vy; 4483 4494 4484 int dummy;4485 4486 4495 /*Get dof list: */ 4487 GetDofList(&doflist [0],&dummy);4496 GetDofList(&doflist); 4488 4497 4489 4498 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ … … 4501 4510 VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); 4502 4511 4512 /*Free ressources:*/ 4513 xfree((void**)&doflist); 4503 4514 } 4504 4515 /*}}}*/ … … 4512 4523 const int numdof=numdofpervertex*numvertices; 4513 4524 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}}; 4514 4515 int doflist[numdof]; 4525 int* doflist=NULL; 4516 4526 double values[numdof]; 4517 4527 double vz; 4518 4528 4519 int dummy;4520 4521 4529 /*Get dof list: */ 4522 GetDofList(&doflist [0],&dummy);4530 GetDofList(&doflist); 4523 4531 4524 4532 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ … … 4533 4541 /*Add value to global vector*/ 4534 4542 VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); 4543 4544 /*Free ressources:*/ 4545 xfree((void**)&doflist); 4535 4546 } 4536 4547 /*}}}*/ … … 4544 4555 const int numdof=numdofpervertex*numvertices; 4545 4556 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}}; 4546 4547 int doflist[numdof]; 4557 int* doflist=NULL; 4548 4558 double values[numdof]; 4549 4559 double vx,vy,vz,p; 4550 4551 int dummy;4552 4560 double stokesreconditioning; 4553 4561 4554 4562 /*Get dof list: */ 4555 GetDofList(&doflist [0],&dummy);4563 GetDofList(&doflist); 4556 4564 4557 4565 /*Recondition pressure: */ … … 4576 4584 VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); 4577 4585 4586 /*Free ressources:*/ 4587 xfree((void**)&doflist); 4578 4588 } 4579 4589 /*}}}*/ … … 4587 4597 const int numdof=numdofpervertex*numvertices; 4588 4598 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}}; 4589 4590 int doflist[numdof]; 4599 int* doflist=NULL; 4591 4600 double values[numdof]; 4592 4601 double vz; 4593 4602 4594 int dummy;4595 4603 4596 4604 /*Get dof list: */ 4597 GetDofList(&doflist [0],&dummy);4605 GetDofList(&doflist); 4598 4606 4599 4607 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ … … 4608 4616 /*Add value to global vector*/ 4609 4617 VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); 4618 4619 /*Free ressources:*/ 4620 xfree((void**)&doflist); 4610 4621 } 4611 4622 /*}}}*/ … … 4907 4918 const int numdofpervertex=2; 4908 4919 const int numdof=numdofpervertex*numvertices; 4909 4910 int doflist[numdof]; 4920 int* doflist=NULL; 4911 4921 double values[numdof]; 4912 4922 double vx[numvertices]; … … 4928 4938 4929 4939 /*Get dof list: */ 4930 GetDofList(&doflist [0],&dummy);4940 GetDofList(&doflist); 4931 4941 4932 4942 /*Use the dof list to index into the solution vector: */ … … 4988 4998 penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id); 4989 4999 } 5000 5001 /*Free ressources:*/ 5002 xfree((void**)&doflist); 4990 5003 } 4991 5004 /*}}}*/ … … 4998 5011 const int numdofpervertex=2; 4999 5012 const int numdof=numdofpervertex*numvertices; 5000 5001 int doflist[numdof]; 5013 int* doflist=NULL; 5002 5014 double values[numdof]; 5003 5015 double vx[numvertices]; … … 5016 5028 5017 5029 /*Get dof list: */ 5018 GetDofList(&doflist [0],&dummy);5030 GetDofList(&doflist); 5019 5031 5020 5032 /*Get node data: */ … … 5068 5080 this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); 5069 5081 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 5082 5083 /*Free ressources:*/ 5084 xfree((void**)&doflist); 5070 5085 } 5071 5086 … … 5079 5094 const int numdofpervertex=2; 5080 5095 const int numdof=numdofpervertex*numvertices; 5081 5082 int doflist[numdof]; 5096 int* doflist=NULL; 5083 5097 double values[numdof]; 5084 5098 double vx[numvertices]; … … 5097 5111 5098 5112 /*Get dof list: */ 5099 GetDofList(&doflist [0],&dummy);5113 GetDofList(&doflist); 5100 5114 5101 5115 /*Get node data: */ … … 5149 5163 this->inputs->AddInput(new TriaVertexInput(VelEnum,vel)); 5150 5164 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 5165 5166 /*Free ressources:*/ 5167 xfree((void**)&doflist); 5151 5168 } 5152 5169 … … 5160 5177 const int numdofpervertex=1; 5161 5178 const int numdof=numdofpervertex*numvertices; 5162 5163 int doflist[numdof]; 5179 int* doflist=NULL; 5164 5180 double values[numdof]; 5165 5181 double vx[numvertices]; … … 5180 5196 5181 5197 /*Get dof list: */ 5182 GetDofList(&doflist [0],&dummy);5198 GetDofList(&doflist); 5183 5199 5184 5200 /*Get node data: */ … … 5233 5249 this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); 5234 5250 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 5235 } /*}}}*/ 5251 5252 /*Free ressources:*/ 5253 xfree((void**)&doflist); 5254 } 5255 /*}}}*/ 5236 5256 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticStokes {{{1*/ 5237 5257 void Penta::InputUpdateFromSolutionDiagnosticStokes(double* solution){ … … 5242 5262 const int numdofpervertex=4; 5243 5263 const int numdof=numdofpervertex*numvertices; 5244 5245 int doflist[numdof]; 5264 int* doflist=NULL; 5246 5265 double values[numdof]; 5247 5266 double vx[numvertices]; … … 5250 5269 double vel[numvertices]; 5251 5270 double pressure[numvertices]; 5252 int dummy;5253 5271 double stokesreconditioning; 5254 5272 5255 5273 /*Get dof list: */ 5256 GetDofList(&doflist [0],&dummy);5274 GetDofList(&doflist); 5257 5275 5258 5276 /*Use the dof list to index into the solution vector: */ … … 5291 5309 this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); 5292 5310 this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); 5311 5312 /*Free ressources:*/ 5313 xfree((void**)&doflist); 5293 5314 } 5294 5315 … … 5302 5323 const int numdofpervertex=4; 5303 5324 const int numdof=numdofpervertex*numvertices; 5304 5305 int doflist[numdof]; 5325 int* doflist=NULL; 5306 5326 double values[numdof]; 5307 5327 double lambdax[numvertices]; … … 5310 5330 double lambdap[numvertices]; 5311 5331 5312 int dummy;5313 5314 5332 /*Get dof list: */ 5315 GetDofList(&doflist [0],&dummy);5333 GetDofList(&doflist); 5316 5334 5317 5335 /*Use the dof list to index into the solution vector: */ … … 5334 5352 this->inputs->AddInput(new PentaVertexInput(AdjointpEnum,lambdap)); 5335 5353 5354 /*Free ressources:*/ 5355 xfree((void**)&doflist); 5336 5356 } 5337 5357 /*}}}*/ … … 5344 5364 const int numdofpervertex=2; 5345 5365 const int numdof=numdofpervertex*numvertices; 5346 5347 int doflist[numdof]; 5366 int* doflist=NULL; 5348 5367 double values[numdof]; 5349 5368 double lambdax[numvertices]; 5350 5369 double lambday[numvertices]; 5351 5370 5352 int dummy;5353 5354 5371 /*Get dof list: */ 5355 GetDofList(&doflist [0],&dummy);5372 GetDofList(&doflist); 5356 5373 5357 5374 /*Use the dof list to index into the solution vector: */ … … 5370 5387 this->inputs->AddInput(new PentaVertexInput(AdjointyEnum,lambday)); 5371 5388 5389 /*Free ressources:*/ 5390 xfree((void**)&doflist); 5372 5391 } 5373 5392 /*}}}*/ … … 5380 5399 const int numdofpervertex=1; 5381 5400 const int numdof=numdofpervertex*numvertices; 5382 5383 int doflist[numdof]; 5401 int* doflist=NULL; 5384 5402 double values[numdof]; 5385 5403 double B[numdof]; 5386 5404 double B_average; 5387 5405 5388 int dummy;5389 5390 5406 /*Get dof list: */ 5391 GetDofList(&doflist [0],&dummy);5407 GetDofList(&doflist); 5392 5408 5393 5409 /*Use the dof list to index into the solution vector: */ … … 5407 5423 this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); 5408 5424 5425 /*Free ressources:*/ 5426 xfree((void**)&doflist); 5409 5427 } 5410 5428 /*}}}*/ … … 5415 5433 const int numdofpervertex = 1; 5416 5434 const int numdof = numdofpervertex *numvertices; 5417 int doflist[numdof];5435 int* doflist=NULL; 5418 5436 double values[numdof]; 5419 int dummy;5420 5437 5421 5438 /*Get dof list: */ 5422 GetDofList(&doflist [0],&dummy);5439 GetDofList(&doflist); 5423 5440 5424 5441 /*Use the dof list to index into the solution vector: */ … … 5429 5446 /*Add input to the element: */ 5430 5447 this->inputs->AddInput(new PentaVertexInput(enum_type,values)); 5448 5449 /*Free ressources:*/ 5450 xfree((void**)&doflist); 5431 5451 } 5432 5452 /*}}}*/ … … 5438 5458 const int numdof = numdofpervertex *numvertices; 5439 5459 const int numdof2d = numdof/2; 5440 int doflist[numdof];5460 int* doflist=NULL; 5441 5461 double values[numdof]; 5442 int dummy;5443 5462 Penta *penta = NULL; 5444 5463 bool onbed; … … 5451 5470 5452 5471 /*Get dof list: */ 5453 GetDofList(&doflist [0],&dummy);5472 GetDofList(&doflist); 5454 5473 5455 5474 /*Use the dof list to index into the solution vector and extrude it */ … … 5472 5491 penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id); 5473 5492 } 5493 5494 /*Free ressources:*/ 5495 xfree((void**)&doflist); 5474 5496 } 5475 5497 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.