![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 15726
- Timestamp:
- 08/06/13 13:44:14 (12 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 4 edited
-
classes/Elements/Penta.cpp (modified) (2 diffs)
-
classes/Elements/Tria.cpp (modified) (2 diffs)
-
classes/FemModel.cpp (modified) (1 diff)
-
modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r15724 r15726 5286 5286 ElementVector* Penta::CreatePVectorAdjointHO(void){ 5287 5287 5288 5289 /*Nothing to be done if not on surface*/ 5288 5290 if (!IsOnSurface()) return NULL; 5289 5291 5290 /*Call Tria function*/ 5291 Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 5292 ElementVector* pe=tria->CreatePVectorAdjointHoriz(); 5293 delete tria->material; delete tria; 5294 5295 /*clean up and return*/ 5292 /*Intermediaries */ 5293 int i,j,resp; 5294 int *responses=NULL; 5295 int num_responses; 5296 IssmDouble Jdet2d; 5297 IssmDouble obs_velocity_mag,velocity_mag; 5298 IssmDouble dux,duy; 5299 IssmDouble epsvel=2.220446049250313e-16; 5300 IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/ 5301 IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.; 5302 IssmDouble vx,vy,vxobs,vyobs,weight; 5303 IssmDouble xyz_list[NUMVERTICES][3]; 5304 IssmDouble xyz_list_tria[NUMVERTICES2D][3]; 5305 5306 /*Fetch number of nodes and dof for this finite element*/ 5307 int numnodes = this->NumberofNodes(); 5308 5309 /*Initialize Element vector*/ 5310 ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); 5311 IssmDouble* basis = xNew<IssmDouble>(numnodes); 5312 5313 /*Retrieve all inputs and parameters*/ 5314 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 5315 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); 5316 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 5317 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 5318 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input); 5319 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input); 5320 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 5321 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 5322 5323 for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j]; 5324 5325 /*Get Surface if required by one response*/ 5326 for(resp=0;resp<num_responses;resp++){ 5327 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 5328 inputs->GetInputValue(&S,SurfaceAreaEnum); break; 5329 } 5330 } 5331 5332 /* Start looping on the number of gaussian points: */ 5333 GaussPenta* gauss=new GaussPenta(3,4,5,4); 5334 for(int ig=gauss->begin();ig<gauss->end();ig++){ 5335 5336 gauss->GaussPoint(ig); 5337 5338 /* Get Jacobian determinant: */ 5339 GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss); 5340 5341 /*Get all parameters at gaussian point*/ 5342 vx_input->GetInputValue(&vx,gauss); 5343 vy_input->GetInputValue(&vy,gauss); 5344 vxobs_input->GetInputValue(&vxobs,gauss); 5345 vyobs_input->GetInputValue(&vyobs,gauss); 5346 GetNodalFunctions(basis, gauss); 5347 5348 /*Loop over all requested responses*/ 5349 for(resp=0;resp<num_responses;resp++){ 5350 5351 weights_input->GetInputValue(&weight,gauss,resp); 5352 5353 switch(responses[resp]){ 5354 case SurfaceAbsVelMisfitEnum: 5355 /* 5356 * 1 [ 2 2 ] 5357 * J = --- | (u - u ) + (v - v ) | 5358 * 2 [ obs obs ] 5359 * 5360 * dJ 5361 * DU = - -- = (u - u ) 5362 * du obs 5363 */ 5364 for(i=0;i<numnodes;i++){ 5365 dux=vxobs-vx; 5366 duy=vyobs-vy; 5367 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; 5368 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; 5369 } 5370 break; 5371 case SurfaceRelVelMisfitEnum: 5372 /* 5373 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 5374 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 5375 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 5376 * obs obs 5377 * 5378 * dJ \bar{v}^2 5379 * DU = - -- = ------------- (u - u ) 5380 * du (u + eps)^2 obs 5381 * obs 5382 */ 5383 for(i=0;i<numnodes;i++){ 5384 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 5385 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 5386 dux=scalex*(vxobs-vx); 5387 duy=scaley*(vyobs-vy); 5388 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; 5389 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; 5390 } 5391 break; 5392 case SurfaceLogVelMisfitEnum: 5393 /* 5394 * [ vel + eps ] 2 5395 * J = 4 \bar{v}^2 | log ( ----------- ) | 5396 * [ vel + eps ] 5397 * obs 5398 * 5399 * dJ 2 * log(...) 5400 * DU = - -- = - 4 \bar{v}^2 ------------- u 5401 * du vel^2 + eps 5402 * 5403 */ 5404 for(i=0;i<numnodes;i++){ 5405 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 5406 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 5407 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 5408 dux=scale*vx; 5409 duy=scale*vy; 5410 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; 5411 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; 5412 } 5413 break; 5414 case SurfaceAverageVelMisfitEnum: 5415 /* 5416 * 1 2 2 5417 * J = --- sqrt( (u - u ) + (v - v ) ) 5418 * S obs obs 5419 * 5420 * dJ 1 1 5421 * DU = - -- = - --- ----------- * 2 (u - u ) 5422 * du S 2 sqrt(...) obs 5423 */ 5424 for(i=0;i<numnodes;i++){ 5425 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 5426 dux=scale*(vxobs-vx); 5427 duy=scale*(vyobs-vy); 5428 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; 5429 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; 5430 } 5431 break; 5432 case SurfaceLogVxVyMisfitEnum: 5433 /* 5434 * 1 [ |u| + eps 2 |v| + eps 2 ] 5435 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 5436 * 2 [ |u |+ eps |v |+ eps ] 5437 * obs obs 5438 * dJ 1 u 1 5439 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 5440 * du |u| + eps |u| u + eps 5441 */ 5442 for(i=0;i<numnodes;i++){ 5443 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 5444 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 5445 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i]; 5446 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i]; 5447 } 5448 break; 5449 case DragCoefficientAbsGradientEnum: 5450 /*Nothing in P vector*/ 5451 break; 5452 case ThicknessAbsGradientEnum: 5453 /*Nothing in P vector*/ 5454 break; 5455 case ThicknessAlongGradientEnum: 5456 /*Nothing in P vector*/ 5457 break; 5458 case ThicknessAcrossGradientEnum: 5459 /*Nothing in P vector*/ 5460 break; 5461 case RheologyBbarAbsGradientEnum: 5462 /*Nothing in P vector*/ 5463 break; 5464 default: 5465 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); 5466 } 5467 } 5468 } 5469 5470 /*Clean up and return*/ 5471 xDelete<IssmDouble>(basis); 5472 xDelete<int>(responses); 5473 delete gauss; 5296 5474 return pe; 5297 5475 } … … 5583 5761 break; 5584 5762 case DragCoefficientAbsGradientEnum: 5585 if (!IsOnBed()) return; 5586 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 5587 tria->GradjDragGradient(gradient,resp,control_index); 5588 delete tria->material; delete tria; 5763 if(IsOnBed()){ 5764 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 5765 tria->GradjDragGradient(gradient,resp,control_index); 5766 delete tria->material; delete tria; 5767 } 5589 5768 break; 5590 5769 case RheologyBbarAbsGradientEnum: 5591 if (!IsOnBed()) return; 5592 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 5593 tria->GradjBGradient(gradient,resp,control_index); 5594 delete tria->material; delete tria; 5770 if(IsOnBed()){ 5771 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 5772 tria->GradjBGradient(gradient,resp,control_index); 5773 delete tria->material; delete tria; 5774 } 5595 5775 break; 5596 5776 default: -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r15718 r15726 6373 6373 D[0][0]=D_scalar*dvxdx; 6374 6374 D[0][1]=0.; 6375 D[1][0]=0.; 6375 6376 D[1][1]=D_scalar*dvydy; 6376 D[1][0]=0.;6377 6377 TripleMultiply(B,2,numnodes,1, 6378 6378 &D[0][0],2,2,0, … … 6823 6823 D[0][0]=D_scalar*dvxdx; 6824 6824 D[0][1]=0.; 6825 D[1][0]=0.; 6825 6826 D[1][1]=D_scalar*dvydy; 6826 D[1][0]=0.;6827 6827 TripleMultiply(B,2,numnodes,1, 6828 6828 &D[0][0],2,2,0, -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r15688 r15726 153 153 /*Finalize PETSC for this model: */ 154 154 #ifdef _HAVE_PETSC_ 155 _printf0_("closing P etsc\n");155 _printf0_("closing PETSc\n"); 156 156 PetscFinalize(); 157 157 #endif -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
r15722 r15726 32 32 iomodel->Constant(&materials_type,MaterialsEnum); 33 33 34 /*Now, is the flag macayaealHO on? otherwise, do nothing: */ 35 if(!isSSA & !isL1L2 & !isHO & !isFS) return; 36 34 37 /*Fetch data needed and allocate vectors: */ 35 38 iomodel->FetchData(1,FlowequationElementEquationEnum); 36 39 finiteelement_list=xNewZeroInit<int>(iomodel->numberofelements); 37 40 38 /*Now, is the flag macayaealHO on? otherwise, do nothing: */39 if(!isSSA & !isL1L2 & !isHO & !isFS) return;40 41 41 42 /*Do we have coupling*/
Note:
See TracChangeset
for help on using the changeset viewer.
