![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 16847
- Timestamp:
- 11/21/13 08:31:55 (11 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 9 edited
-
analyses/AdjointHorizAnalysis.cpp (modified) (1 diff)
-
analyses/AdjointHorizAnalysis.h (modified) (1 diff)
-
classes/Elements/Element.h (modified) (1 diff)
-
classes/Elements/Penta.cpp (modified) (1 diff)
-
classes/Elements/Penta.h (modified) (1 diff)
-
classes/Elements/Seg.cpp (modified) (1 diff)
-
classes/Elements/Seg.h (modified) (1 diff)
-
classes/Elements/Tria.cpp (modified) (1 diff)
-
classes/Elements/Tria.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp ¶
r16782 r16847 30 30 }/*}}}*/ 31 31 ElementVector* AdjointHorizAnalysis::CreatePVector(Element* element){/*{{{*/ 32 _error_("not implemented yet"); 32 33 int approximation; 34 element->GetInputValue(&approximation,ApproximationEnum); 35 switch(approximation){ 36 case SSAApproximationEnum: 37 return CreatePVectorSSA(element); 38 case HOApproximationEnum: 39 return CreatePVectorHO(element); 40 case FSApproximationEnum: 41 return CreatePVectorFS(element); 42 case NoneApproximationEnum: 43 return NULL; 44 default: 45 _error_("Approximation "<<EnumToStringx(approximation)<<" not supported"); 46 } 47 }/*}}}*/ 48 ElementVector* AdjointHorizAnalysis::CreatePVectorFS(Element* element){/*{{{*/ 49 50 /*Nothing to be done if not on surface*/ 51 if(!element->IsOnSurface()) return NULL; 52 53 /*Intermediaries */ 54 int num_responses,i,meshtype,dim; 55 IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag; 56 IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; 57 IssmDouble scalex,scaley,scale,S; 58 int *responses = NULL; 59 IssmDouble *xyz_list_top = NULL; 60 61 /*Get problem dimension*/ 62 element->FindParam(&meshtype,MeshTypeEnum); 63 switch(meshtype){ 64 case Mesh2DverticalEnum: dim = 2; break; 65 case Mesh3DEnum: dim = 3; break; 66 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet"); 67 } 68 69 /*Fetch number of nodes and dof for this finite element*/ 70 int vnumnodes = element->GetNumberOfNodesVelocity(); 71 int pnumnodes = element->GetNumberOfNodesPressure(); 72 73 /*Prepare coordinate system list*/ 74 int* cs_list = xNew<int>(vnumnodes+pnumnodes); 75 for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum; 76 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum; 77 78 /*Initialize Element vector and vectors*/ 79 ElementVector* pe = element->NewElementVector(FSApproximationEnum); 80 IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes); 81 82 /*Retrieve all inputs and parameters*/ 83 element->GetVerticesCoordinatesTop(&xyz_list_top); 84 element->FindParam(&num_responses,InversionNumCostFunctionsEnum); 85 element->FindParam(&responses,NULL,InversionCostFunctionsEnum); 86 Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 87 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 88 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); 89 Input* vxobs_input = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 90 Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 91 IssmDouble epsvel = 2.220446049250313e-16; 92 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ 93 94 /*Get Surface if required by one response*/ 95 for(int resp=0;resp<num_responses;resp++){ 96 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 97 element->GetInputValue(&S,SurfaceAreaEnum); break; 98 } 99 } 100 101 /* Start looping on the number of gaussian points: */ 102 Gauss* gauss=element->NewGaussTop(4); 103 for(int ig=gauss->begin();ig<gauss->end();ig++){ 104 gauss->GaussPoint(ig); 105 106 element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss); 107 element->NodalFunctionsVelocity(vbasis,gauss); 108 109 vx_input->GetInputValue(&vx,gauss); 110 vy_input->GetInputValue(&vy,gauss); 111 vxobs_input->GetInputValue(&vxobs,gauss); 112 vyobs_input->GetInputValue(&vyobs,gauss); 113 114 /*Loop over all requested responses*/ 115 for(int resp=0;resp<num_responses;resp++){ 116 weights_input->GetInputValue(&weight,gauss,responses[resp]); 117 118 switch(responses[resp]){ 119 case SurfaceAbsVelMisfitEnum: 120 /* 121 * 1 [ 2 2 ] 122 * J = --- | (u - u ) + (v - v ) | 123 * 2 [ obs obs ] 124 * 125 * dJ 126 * DU = - -- = (u - u ) 127 * du obs 128 */ 129 for(i=0;i<vnumnodes;i++){ 130 dux=vxobs-vx; 131 duy=vyobs-vy; 132 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 133 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 134 } 135 break; 136 case SurfaceRelVelMisfitEnum: 137 /* 138 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 139 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 140 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 141 * obs obs 142 * 143 * dJ \bar{v}^2 144 * DU = - -- = ------------- (u - u ) 145 * du (u + eps)^2 obs 146 * obs 147 */ 148 for(i=0;i<vnumnodes;i++){ 149 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 150 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 151 dux=scalex*(vxobs-vx); 152 duy=scaley*(vyobs-vy); 153 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 154 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 155 } 156 break; 157 case SurfaceLogVelMisfitEnum: 158 /* 159 * [ vel + eps ] 2 160 * J = 4 \bar{v}^2 | log ( ----------- ) | 161 * [ vel + eps ] 162 * obs 163 * 164 * dJ 2 * log(...) 165 * DU = - -- = - 4 \bar{v}^2 ------------- u 166 * du vel^2 + eps 167 * 168 */ 169 for(i=0;i<vnumnodes;i++){ 170 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 171 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 172 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 173 dux=scale*vx; 174 duy=scale*vy; 175 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 176 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 177 } 178 break; 179 case SurfaceAverageVelMisfitEnum: 180 /* 181 * 1 2 2 182 * J = --- sqrt( (u - u ) + (v - v ) ) 183 * S obs obs 184 * 185 * dJ 1 1 186 * DU = - -- = - --- ----------- * 2 (u - u ) 187 * du S 2 sqrt(...) obs 188 */ 189 for(i=0;i<vnumnodes;i++){ 190 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 191 dux=scale*(vxobs-vx); 192 duy=scale*(vyobs-vy); 193 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 194 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 195 } 196 break; 197 case SurfaceLogVxVyMisfitEnum: 198 /* 199 * 1 [ |u| + eps 2 |v| + eps 2 ] 200 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 201 * 2 [ |u |+ eps |v |+ eps ] 202 * obs obs 203 * dJ 1 u 1 204 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 205 * du |u| + eps |u| u + eps 206 */ 207 for(i=0;i<vnumnodes;i++){ 208 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 209 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 210 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 211 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 212 } 213 break; 214 case DragCoefficientAbsGradientEnum: 215 /*Nothing in P vector*/ 216 break; 217 case ThicknessAbsGradientEnum: 218 /*Nothing in P vector*/ 219 break; 220 case ThicknessAlongGradientEnum: 221 /*Nothing in P vector*/ 222 break; 223 case ThicknessAcrossGradientEnum: 224 /*Nothing in P vector*/ 225 break; 226 case RheologyBbarAbsGradientEnum: 227 /*Nothing in P vector*/ 228 break; 229 default: 230 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); 231 } 232 } 233 } 234 235 /*Transform coordinate system*/ 236 element->TransformLoadVectorCoord(pe,cs_list); 237 238 /*Clean up and return*/ 239 xDelete<int>(cs_list); 240 xDelete<IssmDouble>(xyz_list_top); 241 xDelete<IssmDouble>(vbasis); 242 delete gauss; 243 return pe; 244 _error_("S"); 245 }/*}}}*/ 246 ElementVector* AdjointHorizAnalysis::CreatePVectorHO(Element* element){/*{{{*/ 247 248 /*Nothing to be done if not on surface*/ 249 if(!element->IsOnSurface()) return NULL; 250 251 /*Intermediaries */ 252 int num_responses,i; 253 IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag; 254 IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; 255 IssmDouble scalex,scaley,scale,S; 256 int *responses = NULL; 257 IssmDouble *xyz_list_top = NULL; 258 259 /*Fetch number of nodes and dof for this finite element*/ 260 int numnodes = element->GetNumberOfNodes(); 261 262 /*Initialize Element vector and vectors*/ 263 ElementVector* pe = element->NewElementVector(HOApproximationEnum); 264 IssmDouble* basis = xNew<IssmDouble>(numnodes); 265 266 /*Retrieve all inputs and parameters*/ 267 element->GetVerticesCoordinatesTop(&xyz_list_top); 268 element->FindParam(&num_responses,InversionNumCostFunctionsEnum); 269 element->FindParam(&responses,NULL,InversionCostFunctionsEnum); 270 Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 271 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 272 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); 273 Input* vxobs_input = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 274 Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 275 IssmDouble epsvel = 2.220446049250313e-16; 276 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ 277 278 /*Get Surface if required by one response*/ 279 for(int resp=0;resp<num_responses;resp++){ 280 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 281 element->GetInputValue(&S,SurfaceAreaEnum); break; 282 } 283 } 284 285 /* Start looping on the number of gaussian points: */ 286 Gauss* gauss=element->NewGaussTop(4); 287 for(int ig=gauss->begin();ig<gauss->end();ig++){ 288 gauss->GaussPoint(ig); 289 290 element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss); 291 element->NodalFunctions(basis, gauss); 292 293 vx_input->GetInputValue(&vx,gauss); 294 vy_input->GetInputValue(&vy,gauss); 295 vxobs_input->GetInputValue(&vxobs,gauss); 296 vyobs_input->GetInputValue(&vyobs,gauss); 297 298 /*Loop over all requested responses*/ 299 for(int resp=0;resp<num_responses;resp++){ 300 weights_input->GetInputValue(&weight,gauss,responses[resp]); 301 302 switch(responses[resp]){ 303 case SurfaceAbsVelMisfitEnum: 304 /* 305 * 1 [ 2 2 ] 306 * J = --- | (u - u ) + (v - v ) | 307 * 2 [ obs obs ] 308 * 309 * dJ 310 * DU = - -- = (u - u ) 311 * du obs 312 */ 313 for(i=0;i<numnodes;i++){ 314 dux=vxobs-vx; 315 duy=vyobs-vy; 316 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 317 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 318 } 319 break; 320 case SurfaceRelVelMisfitEnum: 321 /* 322 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 323 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 324 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 325 * obs obs 326 * 327 * dJ \bar{v}^2 328 * DU = - -- = ------------- (u - u ) 329 * du (u + eps)^2 obs 330 * obs 331 */ 332 for(i=0;i<numnodes;i++){ 333 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 334 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 335 dux=scalex*(vxobs-vx); 336 duy=scaley*(vyobs-vy); 337 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 338 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 339 } 340 break; 341 case SurfaceLogVelMisfitEnum: 342 /* 343 * [ vel + eps ] 2 344 * J = 4 \bar{v}^2 | log ( ----------- ) | 345 * [ vel + eps ] 346 * obs 347 * 348 * dJ 2 * log(...) 349 * DU = - -- = - 4 \bar{v}^2 ------------- u 350 * du vel^2 + eps 351 * 352 */ 353 for(i=0;i<numnodes;i++){ 354 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 355 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 356 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 357 dux=scale*vx; 358 duy=scale*vy; 359 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 360 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 361 } 362 break; 363 case SurfaceAverageVelMisfitEnum: 364 /* 365 * 1 2 2 366 * J = --- sqrt( (u - u ) + (v - v ) ) 367 * S obs obs 368 * 369 * dJ 1 1 370 * DU = - -- = - --- ----------- * 2 (u - u ) 371 * du S 2 sqrt(...) obs 372 */ 373 for(i=0;i<numnodes;i++){ 374 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 375 dux=scale*(vxobs-vx); 376 duy=scale*(vyobs-vy); 377 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 378 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 379 } 380 break; 381 case SurfaceLogVxVyMisfitEnum: 382 /* 383 * 1 [ |u| + eps 2 |v| + eps 2 ] 384 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 385 * 2 [ |u |+ eps |v |+ eps ] 386 * obs obs 387 * dJ 1 u 1 388 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 389 * du |u| + eps |u| u + eps 390 */ 391 for(i=0;i<numnodes;i++){ 392 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 393 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 394 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 395 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 396 } 397 break; 398 case DragCoefficientAbsGradientEnum: 399 /*Nothing in P vector*/ 400 break; 401 case ThicknessAbsGradientEnum: 402 /*Nothing in P vector*/ 403 break; 404 case ThicknessAlongGradientEnum: 405 /*Nothing in P vector*/ 406 break; 407 case ThicknessAcrossGradientEnum: 408 /*Nothing in P vector*/ 409 break; 410 case RheologyBbarAbsGradientEnum: 411 /*Nothing in P vector*/ 412 break; 413 default: 414 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); 415 } 416 } 417 } 418 419 /*Transform coordinate system*/ 420 element->TransformLoadVectorCoord(pe,XYEnum); 421 422 /*Clean up and return*/ 423 xDelete<IssmDouble>(xyz_list_top); 424 xDelete<IssmDouble>(basis); 425 delete gauss; 426 return pe; 427 428 }/*}}}*/ 429 ElementVector* AdjointHorizAnalysis::CreatePVectorSSA(Element* element){/*{{{*/ 430 431 /*Intermediaries*/ 432 int meshtype; 433 Element* basalelement; 434 435 /*Get basal element*/ 436 element->FindParam(&meshtype,MeshTypeEnum); 437 switch(meshtype){ 438 case Mesh2DhorizontalEnum: 439 basalelement = element; 440 break; 441 case Mesh3DEnum: 442 if(!element->IsOnBed()) return NULL; 443 basalelement = element->SpawnBasalElement(); 444 break; 445 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet"); 446 } 447 448 /*Intermediaries */ 449 int num_responses,i; 450 IssmDouble thickness,Jdet,obs_velocity_mag,velocity_mag; 451 IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; 452 IssmDouble scalex,scaley,scale,S; 453 int *responses = NULL; 454 IssmDouble *xyz_list = NULL; 455 456 /*Fetch number of nodes and dof for this finite element*/ 457 int numnodes = basalelement->GetNumberOfNodes(); 458 459 /*Initialize Element vector and vectors*/ 460 ElementVector* pe = basalelement->NewElementVector(SSAApproximationEnum); 461 IssmDouble* basis = xNew<IssmDouble>(numnodes); 462 463 /*Retrieve all inputs and parameters*/ 464 basalelement->GetVerticesCoordinates(&xyz_list); 465 basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum); 466 basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum); 467 Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 468 Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); 469 Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); 470 Input* vxobs_input = basalelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 471 Input* vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 472 IssmDouble epsvel = 2.220446049250313e-16; 473 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ 474 475 /*Get Surface if required by one response*/ 476 for(int resp=0;resp<num_responses;resp++){ 477 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 478 basalelement->GetInputValue(&S,SurfaceAreaEnum); break; 479 } 480 } 481 482 /* Start looping on the number of gaussian points: */ 483 Gauss* gauss=basalelement->NewGauss(4); 484 for(int ig=gauss->begin();ig<gauss->end();ig++){ 485 gauss->GaussPoint(ig); 486 487 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 488 basalelement->NodalFunctions(basis, gauss); 489 490 vx_input->GetInputValue(&vx,gauss); 491 vy_input->GetInputValue(&vy,gauss); 492 vxobs_input->GetInputValue(&vxobs,gauss); 493 vyobs_input->GetInputValue(&vyobs,gauss); 494 495 /*Loop over all requested responses*/ 496 for(int resp=0;resp<num_responses;resp++){ 497 weights_input->GetInputValue(&weight,gauss,responses[resp]); 498 499 switch(responses[resp]){ 500 case SurfaceAbsVelMisfitEnum: 501 /* 502 * 1 [ 2 2 ] 503 * J = --- | (u - u ) + (v - v ) | 504 * 2 [ obs obs ] 505 * 506 * dJ 507 * DU = - -- = (u - u ) 508 * du obs 509 */ 510 for(i=0;i<numnodes;i++){ 511 dux=vxobs-vx; 512 duy=vyobs-vy; 513 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 514 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 515 } 516 break; 517 case SurfaceRelVelMisfitEnum: 518 /* 519 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 520 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 521 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 522 * obs obs 523 * 524 * dJ \bar{v}^2 525 * DU = - -- = ------------- (u - u ) 526 * du (u + eps)^2 obs 527 * obs 528 */ 529 for(i=0;i<numnodes;i++){ 530 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 531 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 532 dux=scalex*(vxobs-vx); 533 duy=scaley*(vyobs-vy); 534 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 535 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 536 } 537 break; 538 case SurfaceLogVelMisfitEnum: 539 /* 540 * [ vel + eps ] 2 541 * J = 4 \bar{v}^2 | log ( ----------- ) | 542 * [ vel + eps ] 543 * obs 544 * 545 * dJ 2 * log(...) 546 * DU = - -- = - 4 \bar{v}^2 ------------- u 547 * du vel^2 + eps 548 * 549 */ 550 for(i=0;i<numnodes;i++){ 551 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 552 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 553 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 554 dux=scale*vx; 555 duy=scale*vy; 556 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 557 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 558 } 559 break; 560 case SurfaceAverageVelMisfitEnum: 561 /* 562 * 1 2 2 563 * J = --- sqrt( (u - u ) + (v - v ) ) 564 * S obs obs 565 * 566 * dJ 1 1 567 * DU = - -- = - --- ----------- * 2 (u - u ) 568 * du S 2 sqrt(...) obs 569 */ 570 for(i=0;i<numnodes;i++){ 571 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 572 dux=scale*(vxobs-vx); 573 duy=scale*(vyobs-vy); 574 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 575 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 576 } 577 break; 578 case SurfaceLogVxVyMisfitEnum: 579 /* 580 * 1 [ |u| + eps 2 |v| + eps 2 ] 581 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 582 * 2 [ |u |+ eps |v |+ eps ] 583 * obs obs 584 * dJ 1 u 1 585 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 586 * du |u| + eps |u| u + eps 587 */ 588 for(i=0;i<numnodes;i++){ 589 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 590 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 591 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 592 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 593 } 594 break; 595 case DragCoefficientAbsGradientEnum: 596 /*Nothing in P vector*/ 597 break; 598 case ThicknessAbsGradientEnum: 599 /*Nothing in P vector*/ 600 break; 601 case ThicknessAlongGradientEnum: 602 /*Nothing in P vector*/ 603 break; 604 case ThicknessAcrossGradientEnum: 605 /*Nothing in P vector*/ 606 break; 607 case RheologyBbarAbsGradientEnum: 608 /*Nothing in P vector*/ 609 break; 610 default: 611 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); 612 } 613 } 614 } 615 616 /*Transform coordinate system*/ 617 basalelement->TransformLoadVectorCoord(pe,XYEnum); 618 619 /*Clean up and return*/ 620 xDelete<IssmDouble>(xyz_list); 621 xDelete<IssmDouble>(basis); 622 if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; 623 delete gauss; 624 return pe; 33 625 }/*}}}*/ 34 626 void AdjointHorizAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/ -
TabularUnified issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h ¶
r16782 r16847 23 23 ElementMatrix* CreateKMatrix(Element* element); 24 24 ElementVector* CreatePVector(Element* element); 25 ElementVector* CreatePVectorSSA(Element* element); 26 ElementVector* CreatePVectorHO(Element* element); 27 ElementVector* CreatePVectorFS(Element* element); 25 28 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element); 26 29 void InputUpdateFromSolution(IssmDouble* solution,Element* element); -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Element.h ¶
r16844 r16847 55 55 virtual void FindParam(int* pvalue,int paramenum)=0; 56 56 virtual void FindParam(IssmDouble* pvalue,int paramenum)=0; 57 virtual void FindParam(int** pvalues,int* psize,int paramenum)=0; 57 58 virtual int FiniteElement(void)=0; 58 59 virtual IssmDouble GetMaterialParameter(int enum_in)=0; -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Penta.cpp ¶
r16840 r16847 869 869 void Penta::FindParam(IssmDouble* pvalue,int paramenum){ 870 870 this->parameters->FindParam(pvalue,paramenum); 871 } 872 /*}}}*/ 873 /*FUNCTION Penta::FindParam(int** pvalues,int* psize,int paramenum){{{*/ 874 void Penta::FindParam(int** pvalues,int* psize,int paramenum){ 875 this->parameters->FindParam(pvalues,psize,paramenum); 871 876 } 872 877 /*}}}*/ -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Penta.h ¶
r16840 r16847 78 78 void FindParam(int* pvalue,int paramenum); 79 79 void FindParam(IssmDouble* pvalue,int paramenum); 80 void FindParam(int** pvalues,int* psize,int paramenum); 80 81 int FiniteElement(void); 81 82 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Seg.cpp ¶
r16698 r16847 98 98 } 99 99 /*}}}*/ 100 /*FUNCTION Seg::FindParam(int** pvalues,int* psize,int paramenum){{{*/ 101 void Seg::FindParam(int** pvalues,int* psize,int paramenum){ 102 this->parameters->FindParam(pvalues,psize,paramenum); 103 } 104 /*}}}*/ 100 105 /*FUNCTION Seg::FiniteElement{{{*/ 101 106 int Seg::FiniteElement(void){ -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Seg.h ¶
r16844 r16847 90 90 void FindParam(int* pvalue,int paramenum); 91 91 void FindParam(IssmDouble* pvalue,int paramenum); 92 void FindParam(int** pvalues,int* psize,int paramenum); 92 93 int FiniteElement(void); 93 94 Element* GetBasalElement(void){_error_("not implemented yet");}; -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Tria.cpp ¶
r16846 r16847 848 848 void Tria::FindParam(IssmDouble* pvalue,int paramenum){ 849 849 this->parameters->FindParam(pvalue,paramenum); 850 } 851 /*}}}*/ 852 /*FUNCTION Tria::FindParam(int** pvalues,int* psize,int paramenum){{{*/ 853 void Tria::FindParam(int** pvalues,int* psize,int paramenum){ 854 this->parameters->FindParam(pvalues,psize,paramenum); 850 855 } 851 856 /*}}}*/ -
TabularUnified issm/trunk-jpl/src/c/classes/Elements/Tria.h ¶
r16844 r16847 83 83 void FindParam(int* pvalue,int paramenum); 84 84 void FindParam(IssmDouble* pvalue,int paramenum); 85 void FindParam(int** pvalues,int* psize,int paramenum); 85 86 int FiniteElement(void); 86 87 Element* GetBasalElement(void){_error_("not implemented yet");};
Note:
See TracChangeset
for help on using the changeset viewer.
