Changeset 13721 for issm/trunk-jpl/src/c/classes/FemModel.cpp
- Timestamp:
- 10/17/12 15:59:20 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/FemModel.cpp
r13700 r13721 228 228 } 229 229 /*}}}*/ 230 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/ 231 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){ 232 233 /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use 234 * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several 235 * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the 236 * Slope configuration.*/ 237 238 int found=-1; 239 for(int i=0;i<nummodels;i++){ 240 if (analysis_type_list[i]==configuration_type){ 241 found=i; 242 break; 243 } 244 } 245 if(found!=-1) analysis_counter=found; 246 else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses"); 247 248 /*Now, plug analysis_counter and analysis_type inside the parameters: */ 249 this->parameters->SetParam(analysis_counter,AnalysisCounterEnum); 250 this->parameters->SetParam(analysis_type,AnalysisTypeEnum); 251 this->parameters->SetParam(configuration_type,ConfigurationTypeEnum); 252 253 /*configure elements, loads and nodes, for this new analysis: */ 254 this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters); 255 this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters); 256 this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters); 257 258 #ifdef _HAVE_PETSC_ 259 /*take care of petsc options, that depend on this analysis type (present only after model processor)*/ 260 if(this->parameters->Exist(PetscOptionsStringsEnum)){ 261 PetscOptionsFromAnalysis(this->parameters,analysis_type); 262 if(VerboseSolver()) _pprintLine_(" petsc Options set for analysis type: " << EnumToStringx(analysis_type)); 263 } 264 #endif 265 266 } 267 /*}}}*/ 268 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/ 269 void FemModel::SetCurrentConfiguration(int configuration_type){ 270 this->SetCurrentConfiguration(configuration_type,configuration_type); 271 } 272 /*}}}*/ 230 273 /*FUNCTION FemModel::Solve {{{*/ 231 274 void FemModel::Solve(void){ … … 280 323 /*}}}*/ 281 324 282 /*Numerics: */283 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{*/284 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){285 286 /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use287 * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several288 * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the289 * Slope configuration.*/290 291 int found=-1;292 for(int i=0;i<nummodels;i++){293 if (analysis_type_list[i]==configuration_type){294 found=i;295 break;296 }297 }298 if(found!=-1) analysis_counter=found;299 else _error_("Could not find alias for analysis_type " << EnumToStringx(configuration_type) << " in list of FemModel analyses");300 301 /*Now, plug analysis_counter and analysis_type inside the parameters: */302 this->parameters->SetParam(analysis_counter,AnalysisCounterEnum);303 this->parameters->SetParam(analysis_type,AnalysisTypeEnum);304 this->parameters->SetParam(configuration_type,ConfigurationTypeEnum);305 306 /*configure elements, loads and nodes, for this new analysis: */307 this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);308 this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);309 this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);310 311 #ifdef _HAVE_PETSC_312 /*take care of petsc options, that depend on this analysis type (present only after model processor)*/313 if(this->parameters->Exist(PetscOptionsStringsEnum)){314 PetscOptionsFromAnalysis(this->parameters,analysis_type);315 if(VerboseSolver()) _pprintLine_(" petsc Options set for analysis type: " << EnumToStringx(analysis_type));316 }317 #endif318 319 }320 /*}}}*/321 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{*/322 void FemModel::SetCurrentConfiguration(int configuration_type){323 this->SetCurrentConfiguration(configuration_type,configuration_type);324 }325 /*}}}*/326 327 325 /*Modules:*/ 328 int FemModel::UpdateVertexPositionsx(void){ /*{{{*/326 int FemModel::UpdateVertexPositionsx(void){ /*{{{*/ 329 327 330 328 int i; … … 378 376 NodesDofx(nodes,parameters,analysis_type); 379 377 380 }381 /*}}}*/382 void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/383 384 int i;385 386 /*output: */387 IssmDouble dt;388 389 /*intermediary: */390 Element *element = NULL;391 IssmDouble min_dt = 0;392 IssmDouble node_min_dt = 0;393 394 /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */395 element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt();396 397 for (i=1;i<elements->Size();i++){398 element=(Element*)elements->GetObjectByOffset(i);399 dt=element->TimeAdapt();400 if(dt<min_dt)min_dt=dt;401 }402 403 /*Figure out minimum across the cluster: */404 #ifdef _HAVE_MPI_405 MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() );406 MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm());407 min_dt=node_min_dt;408 #endif409 410 /*Assign output pointers:*/411 *pdt=min_dt;412 378 } 413 379 /*}}}*/ … … 440 406 case MaxVzEnum: MaxVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 441 407 case MaxAbsVzEnum: MaxAbsVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 442 case MassFluxEnum: MassFluxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;408 case MassFluxEnum: this->MassFluxx( responses,process_units); break; 443 409 case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 444 410 case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; … … 543 509 } 544 510 /*}}}*/ 511 void FemModel::TimeAdaptx(IssmDouble* pdt){/*{{{*/ 512 513 int i; 514 515 /*output: */ 516 IssmDouble dt; 517 518 /*intermediary: */ 519 Element *element = NULL; 520 IssmDouble min_dt = 0; 521 IssmDouble node_min_dt = 0; 522 523 /*Go through elements, and figure out the minimum of the time steps for each element (using CFL criterion): */ 524 element=(Element*)elements->GetObjectByOffset(0); min_dt=element->TimeAdapt(); 525 526 for (i=1;i<elements->Size();i++){ 527 element=(Element*)elements->GetObjectByOffset(i); 528 dt=element->TimeAdapt(); 529 if(dt<min_dt)min_dt=dt; 530 } 531 532 /*Figure out minimum across the cluster: */ 533 #ifdef _HAVE_MPI_ 534 MPI_Reduce (&min_dt,&node_min_dt,1,MPI_DOUBLE,MPI_MIN,0,IssmComm::GetComm() ); 535 MPI_Bcast(&node_min_dt,1,MPI_DOUBLE,0,IssmComm::GetComm()); 536 min_dt=node_min_dt; 537 #endif 538 539 /*Assign output pointers:*/ 540 *pdt=min_dt; 541 } 542 /*}}}*/ 543 #ifdef _HAVE_RESPONSES_ 544 void FemModel::MassFluxx(IssmDouble* pmass_flux,bool process_units){/*{{{*/ 545 546 int i,j; 547 Element *element = NULL; 548 int element_id; 549 bool ispresent = false; 550 IssmDouble mass_flux = 0; 551 IssmDouble all_mass_flux = 0; 552 int counter; 553 IssmDouble **array = NULL; 554 int M; 555 int *mdims_array = NULL; 556 int *ndims_array = NULL; 557 IssmDouble *segments = NULL; 558 int num_segments; 559 560 /*First, figure out which segment to compute our mass flux on. Start with retrieving qmu_mass_flux_segments: */ 561 this->parameters->FindParam(&ispresent,MassFluxSegmentsPresentEnum); 562 if(!ispresent)_error_("no mass flux segments available!"); 563 this->parameters->FindParam(&array,&M,&mdims_array,&ndims_array,MassFluxSegmentsEnum); 564 565 /*Retrieve index of segments being used for MassFlux computation: */ 566 parameters->FindParam(&counter,IndexEnum); 567 568 /*retrieve segments from array: */ 569 segments = array[counter-1]; //matlab to "C" indexing 570 num_segments = mdims_array[counter-1]; 571 572 /*Go through segments, and then elements, and figure out which elements belong to a segment. 573 * When we find one, use the element to compute the mass flux on the segment: */ 574 for(i=0;i<num_segments;i++){ 575 element_id=reCast<int,IssmDouble>(*(segments+5*i+4)); 576 for(j=0;j<elements->Size();j++){ 577 element=(Element*)this->elements->GetObjectByOffset(j); 578 if (element->Id()==element_id){ 579 /*We found the element which owns this segment, use it to compute the mass flux: */ 580 mass_flux+=element->MassFlux(segments+5*i+0,process_units); 581 break; 582 } 583 } 584 } 585 586 #ifdef _HAVE_MPI_ 587 MPI_Allreduce ( (void*)&mass_flux,(void*)&all_mass_flux,1,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm()); 588 mass_flux=all_mass_flux; 589 #endif 590 591 /*Free ressources:*/ 592 for(i=0;j<M;i++){ 593 IssmDouble* matrix=array[j]; 594 xDelete<IssmDouble>(matrix); 595 } 596 xDelete<int>(mdims_array); 597 xDelete<int>(ndims_array); 598 xDelete<IssmDouble*>(array); 599 600 /*Assign output pointers: */ 601 *pmass_flux=mass_flux; 602 603 }/*}}}*/ 604 #endif 545 605 #ifdef _HAVE_CONTROL_ 546 606 void FemModel::ThicknessAbsGradientx( IssmDouble* pJ, bool process_units, int weight_index){/*{{{*/ … … 601 661 /*}}}*/ 602 662 #endif 603 604 663 #ifdef _HAVE_DAKOTA_ 605 664 void FemModel::DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses){/*{{{*/
Note:
See TracChangeset
for help on using the changeset viewer.