Changeset 3620
- Timestamp:
- 04/26/10 16:46:08 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 6 added
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Beam.cpp
r3612 r3620 11 11 #include "stdio.h" 12 12 #include "./Beam.h" 13 #include "./BeamVertexInput.h" 13 14 #include <string.h> 14 15 #include "../EnumDefinitions/EnumDefinitions.h" … … 21 22 /*FUNCTION Beam::Beam(){{{1*/ 22 23 Beam::Beam(){ 24 this->inputs=NULL; 23 25 return; 24 26 } 25 27 /*}}}*/ 26 28 /*FUNCTION Beam::Beam(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id, ElementProperties* properties){{{1*/ 27 Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id , ElementProperties* beamproperties):29 Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id): 28 30 hnodes(beam_node_ids,2), 29 31 hmatice(&beam_matice_id,1), 30 32 hmatpar(&beam_matpar_id,1), 31 hnumpar(&beam_numpar_id,1), 32 properties(beamproperties) 33 hnumpar(&beam_numpar_id,1) 33 34 { 34 35 35 36 /*all the initialization has been done by the initializer, just fill in the id: */ 36 37 this->id=beam_id; 37 38 return; 38 this->inputs=new Inputs(); 39 39 } 40 40 /*}}}*/ 41 41 /*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, ElementProperties* properties){{{1*/ 42 Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, ElementProperties* beam_properties):42 Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs): 43 43 hnodes(beam_hnodes), 44 44 hmatice(beam_hmatice), 45 45 hmatpar(beam_hmatpar), 46 hnumpar(beam_hnumpar), 47 properties(beam_properties) 46 hnumpar(beam_hnumpar) 48 47 { 49 48 50 49 /*all the initialization has been done by the initializer, just fill in the id: */ 51 50 this->id=beam_id; 52 51 if(beam_inputs){ 52 this->inputs=(Inputs*)beam_inputs->Copy(); 53 } 54 else{ 55 this->inputs=new Inputs(); 56 } 53 57 return; 54 58 } 55 59 /*}}}*/ 56 60 /*FUNCTION Beam::Beam(int i, IoModel* iomodel){{{1*/ 57 Beam::Beam(int i,IoModel* iomodel){ 58 61 Beam::Beam(int index,IoModel* iomodel){ 62 63 64 int i; 59 65 60 66 /*beam constructor input: */ … … 64 70 int beam_numpar_id; 65 71 int beam_node_ids[2]; 66 double beam_h[2]; 67 double beam_s[2]; 68 double beam_b[2]; 69 double beam_k[2]; 70 int beam_onbed; 71 72 72 double nodeinputs[2]; 73 73 74 /*id: */ 74 beam_id=i +1;75 beam_id=index+1; 75 76 this->id=beam_id; 76 77 77 78 /*hooks: */ 78 beam_matice_id=i +1; //refers to the corresponding material property card79 beam_matice_id=index+1; //refers to the corresponding material property card 79 80 beam_matpar_id=iomodel->numberofvertices2d*(iomodel->numlayers-1)+1;//refers to the corresponding matpar property card 80 81 beam_numpar_id=1; 81 beam_node_ids[0]=i +1;82 beam_node_ids[1]=(int)iomodel->uppernodes[i ]; //grid that lays right on top82 beam_node_ids[0]=index+1; 83 beam_node_ids[1]=(int)iomodel->uppernodes[index]; //grid that lays right on top 83 84 84 85 this->hnodes.Init(beam_node_ids,2); … … 87 88 this->hnumpar.Init(&beam_numpar_id,1); 88 89 89 /*properties: */ 90 beam_h[0]=iomodel->thickness[i]; 91 beam_h[1]=iomodel->thickness[(int)(iomodel->uppernodes[i]-1)]; 92 beam_s[0]=iomodel->surface[i]; 93 beam_s[1]=iomodel->surface[(int)(iomodel->uppernodes[i]-1)]; 94 beam_b[0]=iomodel->bed[i]; 95 beam_b[1]=iomodel->bed[(int)(iomodel->uppernodes[i]-1)]; 96 beam_k[0]=iomodel->drag[i]; 97 beam_k[1]=iomodel->drag[(int)(iomodel->uppernodes[i]-1)]; 98 beam_onbed=(bool)iomodel->gridonbed[i]; 99 100 this->properties.Init(2,beam_h, beam_s, beam_b, beam_k,NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, beam_onbed, UNDEF, UNDEF, UNDEF, UNDEF); 90 //intialize inputs, and add as many inputs per element as requested: 91 this->inputs=new Inputs(); 92 93 if (iomodel->thickness) { 94 nodeinputs[0]=iomodel->thickness[index]; 95 nodeinputs[1]=iomodel->thickness[(int)(iomodel->uppernodes[index]-1)]; 96 this->inputs->AddInput(new BeamVertexInput(ThicknessEnum,nodeinputs)); 97 } 98 if (iomodel->surface) { 99 nodeinputs[0]=iomodel->surface[index]; 100 nodeinputs[1]=iomodel->surface[(int)(iomodel->uppernodes[index]-1)]; 101 this->inputs->AddInput(new BeamVertexInput(SurfaceEnum,nodeinputs)); 102 } 103 if (iomodel->bed) { 104 nodeinputs[0]=iomodel->bed[index]; 105 nodeinputs[1]=iomodel->bed[(int)(iomodel->uppernodes[index]-1)]; 106 this->inputs->AddInput(new BeamVertexInput(BedEnum,nodeinputs)); 107 } 108 if (iomodel->drag_coefficient) { 109 nodeinputs[0]=iomodel->drag_coefficient[index]; 110 nodeinputs[1]=iomodel->drag_coefficient[(int)(iomodel->uppernodes[index]-1)]; 111 this->inputs->AddInput(new BeamVertexInput(DragCoefficientEnum,nodeinputs)); 112 } 113 if (iomodel->gridonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->gridonbed[index])); 114 101 115 } 102 116 /*}}}*/ 103 117 /*FUNCTION Beam::~Beam(){{{1*/ 104 118 Beam::~Beam(){ 119 delete inputs; 105 120 return; 106 121 } … … 136 151 Object* Beam::copy() { 137 152 138 return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar, &this->properties);153 return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); 139 154 140 155 } … … 149 164 hmatpar.DeepEcho(); 150 165 hnumpar.DeepEcho(); 151 properties.DeepEcho(); 166 printf(" inputs\n"); 167 inputs->DeepEcho(); 152 168 153 169 return; … … 173 189 hnumpar.Demarshall(&marshalled_dataset); 174 190 175 /*demarshall properties: */176 properties.Demarshall(&marshalled_dataset);191 /*demarshall inputs: */ 192 inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 177 193 178 194 /*return: */ … … 190 206 hmatpar.Echo(); 191 207 hnumpar.Echo(); 192 properties.Echo(); 208 printf(" inputs\n"); 209 inputs->Echo(); 193 210 194 211 return; … … 200 217 char* marshalled_dataset=NULL; 201 218 int enum_type=0; 219 char* marshalled_inputs=NULL; 220 int marshalled_inputs_size; 202 221 203 222 /*recover marshalled_dataset: */ … … 219 238 hnumpar.Marshall(&marshalled_dataset); 220 239 221 /*Marshall properties: */ 222 properties.Marshall(&marshalled_dataset); 240 /*Marshall inputs: */ 241 marshalled_inputs_size=inputs->MarshallSize(); 242 marshalled_inputs=inputs->Marshall(); 243 memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); 244 marshalled_dataset+=marshalled_inputs_size; 245 246 xfree((void**)&marshalled_inputs); 223 247 224 248 *pmarshalled_dataset=marshalled_dataset; … … 234 258 +hmatpar.MarshallSize() 235 259 +hnumpar.MarshallSize() 236 + properties.MarshallSize()260 +inputs->MarshallSize() 237 261 +sizeof(int); //sizeof(int) for enum type 238 262 } … … 259 283 int doflist[numgrids]; 260 284 double pressure[numgrids]; 285 double surface[numgrids]; 261 286 double rho_ice,g; 262 287 double xyz_list[numgrids][3]; 288 double gauss[numgrids][numgrids]={{1,0},{0,1}}; 263 289 264 290 /*dynamic objects pointed to by hooks: */ … … 286 312 rho_ice=matpar->GetRhoIce(); 287 313 g=matpar->GetG(); 314 315 /*recover value of surface at gauss points (0,1) and (1,0): */ 316 inputs->GetParameterValues(&surface[0],&gauss[0][0],2,SurfaceEnum); 288 317 for(i=0;i<numgrids;i++){ 289 pressure[i]=rho_ice*g*( this->properties.s[i]-xyz_list[i][2]);290 } 291 318 pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); 319 } 320 292 321 /*plug local pressure values into global pressure vector: */ 293 322 VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES); … … 303 332 /*}}}*/ 304 333 /*FUNCTION Beam::CostFunction{{{1*/ 305 double Beam::CostFunction( void*,int,int){334 double Beam::CostFunction(int,int){ 306 335 ISSMERROR(" not supported yet!"); 307 336 } … … 316 345 if (sub_analysis_type==HutterAnalysisEnum) { 317 346 318 CreateKMatrixDiagnosticHutter( Kgg, inputs,analysis_type,sub_analysis_type);347 CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type); 319 348 } 320 349 else … … 338 367 double Ke_gg[numdofs][numdofs]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}; 339 368 int numberofdofspernode; 340 369 bool onbed; 370 371 372 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 373 341 374 GetDofList(&doflist[0],&numberofdofspernode); 342 375 343 if ( this->properties.onbed){376 if (onbed){ 344 377 Ke_gg[0][0]=1; 345 378 Ke_gg[1][1]=1; … … 367 400 if (analysis_type==DiagnosticAnalysisEnum) { 368 401 if (sub_analysis_type==HutterAnalysisEnum) { 369 CreatePVectorDiagnosticHutter( pg, inputs,analysis_type,sub_analysis_type);402 CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type); 370 403 } 371 404 else … … 403 436 double* gauss_weights=NULL; 404 437 double gauss_weight; 438 double gauss1[2]={1,0}; 405 439 int ig; 406 440 double l1l2[2]; … … 412 446 double Jdet; 413 447 double ub,vb; 448 double surface,thickness; 449 450 bool onbed; 414 451 415 452 /*dynamic objects pointed to by hooks: */ … … 419 456 Numpar* numpar=NULL; 420 457 421 ParameterInputs* inputs=NULL;422 423 /*recover pointers: */424 inputs=(ParameterInputs*)vinputs;425 426 458 /*recover objects from hooks: */ 427 459 nodes=(Node**)hnodes.deliverp(); … … 433 465 GetDofList(&doflist[0],&numberofdofspernode); 434 466 467 /*recover some inputs: */ 468 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 469 435 470 /*recover parameters: */ 436 471 rho_ice=matpar->GetRhoIce(); … … 440 475 441 476 //recover extra inputs 442 found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,1,(void**)&nodes[0]); //only recover from fist node, second node will have the same slope 443 if(!found)ISSMERROR(" surfaceslopex missing from inputs!"); 444 445 found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,1,(void**)&nodes[0]);//only recover from fist node, second node will have the same slope 446 if(!found)ISSMERROR(" surfaceslopey missing from inputs!"); 477 inputs->GetParameterValue(&slope[0],&gauss1[0],SurfaceSlopexEnum); 478 inputs->GetParameterValue(&slope[1],&gauss1[0],SurfaceSlopeyEnum); 479 inputs->GetParameterValue(&surface,&gauss1[0],SurfaceEnum); 480 inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum); 447 481 448 482 //Get all element grid data: … … 474 508 475 509 for(j=0;j<NDOF2;j++){ 476 pe_g_gaussian[NDOF2+j]=constant_part*pow(( this->properties.s[0]-z_g)/B,n)*slope[j]*Jdet*gauss_weight;510 pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight; 477 511 } 478 512 … … 484 518 485 519 //Deal with lower surface 486 if ( this->properties.onbed){520 if (onbed){ 487 521 488 522 //compute ub 489 constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thi s->properties.h[0];523 constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness; 490 524 ub=constant_part*slope[0]; 491 525 vb=constant_part*slope[1]; … … 505 539 /*}}}*/ 506 540 /*FUNCTION Beam::Du{{{1*/ 507 void Beam::Du( _p_Vec*,void*,int,int){541 void Beam::Du(Vec,int,int){ 508 542 ISSMERROR(" not supported yet!"); 509 543 } … … 654 688 /*}}}*/ 655 689 /*FUNCTION Beam::Gradj{{{1*/ 656 void Beam::Gradj( _p_Vec*, void*, int, int,char*){690 void Beam::Gradj(Vec, int, int,char*){ 657 691 ISSMERROR(" not supported yet!"); 658 692 } 659 693 /*}}}*/ 660 694 /*FUNCTION Beam::GradjB{{{1*/ 661 void Beam::GradjB( _p_Vec*, void*, int, int){695 void Beam::GradjB(Vec, int, int){ 662 696 ISSMERROR(" not supported yet!"); 663 697 } 664 698 /*}}}*/ 665 699 /*FUNCTION Beam::GradjDrag{{{1*/ 666 void Beam::GradjDrag( _p_Vec*, void*, int,int ){700 void Beam::GradjDrag(Vec, int,int ){ 667 701 ISSMERROR(" not supported yet!"); 668 702 } … … 674 708 /*}}}*/ 675 709 /*FUNCTION Beam::Misfit{{{1*/ 676 double Beam::Misfit( void*,int,int){710 double Beam::Misfit(int,int){ 677 711 ISSMERROR(" not supported yet!"); 678 712 } … … 685 719 /*}}}*/ 686 720 /*FUNCTION Beam::SurfaceArea{{{1*/ 687 double Beam::SurfaceArea( void*,int,int){721 double Beam::SurfaceArea(int,int){ 688 722 ISSMERROR(" not supported yet!"); 689 723 } -
issm/trunk/src/c/objects/Beam.h
r3612 r3620 12 12 #include "./Matpar.h" 13 13 #include "../ModelProcessorx/IoModel.h" 14 #include "../DataSet/Inputs.h" 14 15 #include "./Hook.h" 15 16 … … 30 31 Hook hmatpar; //hook to 1 matpar 31 32 Hook hnumpar; //hook to 1 numpar 32 33 33 Inputs* inputs; 34 34 … … 52 52 int GetId(); 53 53 int MyRank(); 54 void Configure( void* loads,void* nodes,void* materials,void* parameters);54 void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); 55 55 Object* copy(); 56 56 void SetClone(int* minranks); … … 78 78 void GetBedList(double*); 79 79 void GetThicknessList(double* thickness_list); 80 void Du( _p_Vec*,void*, int,int);81 void Gradj( _p_Vec*, void*,int, int,char*);82 void GradjDrag( _p_Vec*, void*,int,int );83 void GradjB( _p_Vec*, void*,int,int );84 double Misfit( void*,int,int);85 double SurfaceArea( void*,int,int);86 double CostFunction( void*,int,int);80 void Du(Vec, int,int); 81 void Gradj(Vec, int, int,char*); 82 void GradjDrag(Vec, int,int ); 83 void GradjB(Vec, int,int ); 84 double Misfit(int,int); 85 double SurfaceArea(int,int); 86 double CostFunction(int,int); 87 87 void GetNodalFunctions(double* l1l2, double gauss_coord); 88 88 void GetParameterValue(double* pvalue, double* value_list,double gauss_coord); -
issm/trunk/src/c/objects/Element.h
r3612 r3620 10 10 11 11 #include "./Object.h" 12 #include "../DataSet/DataSet.h" 13 #include "../DataSet/Parameters.h" 12 14 #include "../toolkits/toolkits.h" 13 15 … … 18 20 virtual ~Element(){}; 19 21 virtual int Enum()=0; 20 virtual void Configure( void* loads,void* nodes,void* materials,void* parameters)=0;22 virtual void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters)=0; 21 23 22 24 virtual void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type)=0; -
issm/trunk/src/c/objects/FemModel.cpp
r3496 r3620 39 39 /*FUNCTION FemModel::FemModel {{{1*/ 40 40 FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_vertices, DataSet* femmodel_constraints,DataSet* femmodel_loads, 41 DataSet* femmodel_materials, DataSet* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg,41 DataSet* femmodel_materials,Parameters* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg, 42 42 Mat femmodel_Rmg,Mat femmodel_Gmn,NodeSets* femmodel_nodesets,Vec femmodel_ys,Vec femmodel_ys0){ 43 43 … … 266 266 /*}}}*/ 267 267 /*FUNCTION FemModel::get_parameters {{{1*/ 268 DataSet* FemModel::get_parameters(void){268 Parameters* FemModel::get_parameters(void){ 269 269 return parameters; 270 270 } -
issm/trunk/src/c/objects/FemModel.h
r3463 r3620 11 11 #include "./Object.h" 12 12 #include "../DataSet/DataSet.h" 13 #include "../DataSet/Parameters.h" 13 14 #include "./DofVec.h" 14 15 #include "../toolkits/toolkits.h" … … 27 28 DataSet* loads; 28 29 DataSet* materials; 29 DataSet* parameters;30 Parameters* parameters; 30 31 31 32 DofVec* partition; … … 41 42 FemModel(); 42 43 ~FemModel(); 43 FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials, DataSet* parameters,44 FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,Parameters* parameters, 44 45 DofVec* partition,DofVec* tpartition,DofVec* yg,Mat Rmg,Mat Gmn,NodeSets* nodesets,Vec ys,Vec ys0); 45 46 … … 69 70 DataSet* get_loads(void); 70 71 DataSet* get_materials(void); 71 DataSet* get_parameters(void);72 Parameters* get_parameters(void); 72 73 DofVec* get_partition(void); 73 74 DofVec* get_tpartition(void); -
issm/trunk/src/c/objects/Model.cpp
r3567 r3620 47 47 DataSet* loads=NULL; 48 48 DataSet* materials=NULL; 49 DataSet* parameters=NULL;49 Parameters* parameters=NULL; 50 50 DofVec* partition=NULL; 51 51 DofVec* tpartition=NULL; … … 121 121 DataSet* loads=NULL; 122 122 DataSet* materials=NULL; 123 DataSet* parameters=NULL;123 Parameters* parameters=NULL; 124 124 DofVec* partition=NULL; 125 125 DofVec* tpartition=NULL; -
issm/trunk/src/c/objects/Penta.cpp
r3613 r3620 10 10 11 11 #include "stdio.h" 12 #include "./PentaVertexInput.h" 12 13 #include "./Penta.h" 13 14 #include <string.h> … … 20 21 /*FUNCTION Penta default constructor {{{1*/ 21 22 Penta::Penta(){ 23 this->inputs=NULL; 22 24 return; 23 25 } 24 26 /*}}}*/ 25 27 /*FUNCTION Penta constructor {{{1*/ 26 Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id , ElementProperties* pentaproperties):28 Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id): 27 29 hnodes(penta_node_ids,6), 28 30 hmatice(&penta_matice_id,1), 29 31 hmatpar(&penta_matpar_id,1), 30 hnumpar(&penta_numpar_id,1), 31 properties(pentaproperties) 32 hnumpar(&penta_numpar_id,1) 32 33 { 33 34 34 35 /*all the initialization has been done by the initializer, just fill in the id: */ 35 36 this->id=penta_id; 36 37 return; 37 this->inputs=new Inputs(); 38 38 39 } 39 40 /*}}}*/ 40 41 /*FUNCTION Penta other constructor {{{1*/ 41 Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties):42 Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* penta_inputs): 42 43 hnodes(penta_hnodes), 43 44 hmatice(penta_hmatice), 44 45 hmatpar(penta_hmatpar), 45 hnumpar(penta_hnumpar), 46 properties(penta_properties) 46 hnumpar(penta_hnumpar) 47 47 { 48 48 49 49 /*all the initialization has been done by the initializer, just fill in the id: */ 50 50 this->id=penta_id; 51 52 return; 51 if(penta_inputs){ 52 this->inputs=(Inputs*)penta_inputs->Copy(); 53 } 54 else{ 55 this->inputs=new Inputs(); 56 } 53 57 } 54 58 /*}}}*/ 55 59 /*FUNCTION Penta other constructor {{{1*/ 56 Penta::Penta(int i , IoModel* iomodel){ //i is the element index57 58 int j;60 Penta::Penta(int index, IoModel* iomodel){ //i is the element index 61 62 IssmInt i; 59 63 int penta_node_ids[6]; 60 64 int penta_matice_id; 61 65 int penta_matpar_id; 62 66 int penta_numpar_id; 67 double nodeinputs[6]; 68 bool collapse; 63 69 64 70 /*id: */ 65 this->id=i +1;71 this->id=index+1; 66 72 67 73 /*hooks: */ 68 for( j=0;j<6;j++){ //go recover node ids, needed to initialize the node hook.69 penta_node_ids[ j]=(int)*(iomodel->elements+6*i+j); //ids for vertices are in the elements array from Matlab70 } 71 penta_matice_id=i +1; //refers to the corresponding ice material object74 for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. 75 penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab 76 } 77 penta_matice_id=index+1; //refers to the corresponding ice material object 72 78 penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object 73 79 penta_numpar_id=1; //refers to numerical parameters object 74 80 75 this->hnodes.Init( penta_node_ids,6);81 this->hnodes.Init(&penta_node_ids[0],6); 76 82 this->hmatice.Init(&penta_matice_id,1); 77 83 this->hmatpar.Init(&penta_matpar_id,1); 78 84 this->hnumpar.Init(&penta_numpar_id,1); 79 this->properties.Init(6,penta_node_ids,this->id,iomodel); 85 86 //intialize inputs, and add as many inputs per element as requested: 87 this->inputs=new Inputs(); 88 89 if (iomodel->thickness) { 90 for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1]; 91 this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); 92 } 93 if (iomodel->surface) { 94 for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_node_ids[i]-1]; 95 this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); 96 } 97 if (iomodel->bed) { 98 for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_node_ids[i]-1]; 99 this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs)); 100 } 101 if (iomodel->drag_coefficient) { 102 for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_node_ids[i]-1]; 103 this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs)); 104 105 if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index])); 106 if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); 107 this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); 108 109 } 110 if (iomodel->melting_rate) { 111 for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_node_ids[i]-1]; 112 this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); 113 } 114 if (iomodel->accumulation_rate) { 115 for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_node_ids[i]-1]; 116 this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); 117 } 118 if (iomodel->geothermalflux) { 119 for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_node_ids[i]-1]; 120 this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs)); 121 } 122 123 if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); 124 if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); 125 if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); 126 if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); 80 127 81 128 //elements of type 3 are macayeal type penta. we collapse the formulation on their base. 82 129 if(iomodel->elements_type){ 83 130 if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ 84 this->properties.collapse=1;131 collapse=1; 85 132 } 86 133 else{ 87 this->properties.collapse=0;134 collapse=0; 88 135 } 89 136 } 137 this->inputs->AddInput(new BoolInput(CollapseEnum,(IssmBool)collapse)); 90 138 91 139 } … … 93 141 /*FUNCTION Penta destructor {{{1*/ 94 142 Penta::~Penta(){ 143 delete inputs; 95 144 return; 96 145 } … … 125 174 /*FUNCTION copy {{{1*/ 126 175 Object* Penta::copy() { 127 return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar, &this->properties);176 return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); 128 177 } 129 178 /*}}}*/ … … 148 197 hnumpar.Demarshall(&marshalled_dataset); 149 198 150 /*demarshall properties: */151 properties.Demarshall(&marshalled_dataset);199 /*demarshall inputs: */ 200 inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 152 201 153 202 /*return: */ … … 166 215 hmatpar.DeepEcho(); 167 216 hnumpar.DeepEcho(); 168 properties.DeepEcho(); 217 printf(" inputs\n"); 218 inputs->DeepEcho(); 169 219 170 220 return; … … 181 231 hmatpar.Echo(); 182 232 hnumpar.Echo(); 183 properties.Echo(); 184 185 return; 233 printf(" inputs\n"); 234 inputs->Echo(); 235 236 } 237 /*}}}*/ 238 /*FUNCTION Enum {{{1*/ 239 int Penta::Enum(void){ 240 241 return PentaEnum; 242 186 243 } 187 244 /*}}}*/ … … 191 248 char* marshalled_dataset=NULL; 192 249 int enum_type=0; 250 char* marshalled_inputs=NULL; 251 int marshalled_inputs_size; 193 252 194 253 /*recover marshalled_dataset: */ … … 210 269 hnumpar.Marshall(&marshalled_dataset); 211 270 212 /*Marshall properties: */ 213 properties.Marshall(&marshalled_dataset); 271 /*Marshall inputs: */ 272 marshalled_inputs_size=inputs->MarshallSize(); 273 marshalled_inputs=inputs->Marshall(); 274 memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); 275 marshalled_dataset+=marshalled_inputs_size; 276 277 xfree((void**)&marshalled_inputs); 278 214 279 215 280 *pmarshalled_dataset=marshalled_dataset; … … 225 290 +hmatpar.MarshallSize() 226 291 +hnumpar.MarshallSize() 227 + properties.MarshallSize()292 +inputs->MarshallSize() 228 293 +sizeof(int); //sizeof(int) for enum type 229 294 } … … 240 305 Hook* tria_hmatpar=NULL; 241 306 Hook* tria_hnumpar=NULL; 242 ElementProperties* tria_properties=NULL; 243 DataSet* tria_inputs=NULL; 307 Inputs* tria_inputs=NULL; 244 308 245 309 indices[0]=g0; … … 251 315 tria_hmatpar=this->hmatpar.Spawn(&zero,1); 252 316 tria_hnumpar=this->hnumpar.Spawn(&zero,1); 253 tria_properties=(ElementProperties*)this->properties.Spawn(indices,3); 254 tria_inputs=(DataSet*)this->inputs->Spawn(indices,3); 255 256 tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_properties,tria_inputs); 317 tria_inputs=(Inputs*)this->inputs->Spawn(indices,3); 318 319 tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_inputs); 257 320 258 321 delete tria_hnodes; … … 260 323 delete tria_hmatpar; 261 324 delete tria_hnumpar; 262 delete tria_properties;263 325 delete tria_inputs; 264 326 … … 266 328 } 267 329 /*}}}*/ 330 331 /*updates: */ 268 332 /*FUNCTION UpdateFromDakota {{{1*/ 269 333 void Penta::UpdateFromDakota(void* vinputs){ … … 275 339 /*FUNCTION Penta::UpdateInputs {{{1*/ 276 340 void Penta::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){ 341 342 /*Just branch to the correct UpdateInputs generator, according to the type of analysis we are carrying out: */ 343 if (analysis_type==ControlAnalysisEnum){ 344 345 UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type); 346 } 347 else if (analysis_type==DiagnosticAnalysisEnum){ 348 349 if (sub_analysis_type==HorizAnalysisEnum){ 350 351 UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type); 352 } 353 else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); 354 355 } 356 else if (analysis_type==SlopecomputeAnalysisEnum){ 357 358 UpdateInputsSlopeCompute( solution,analysis_type,sub_analysis_type); 359 } 360 else if (analysis_type==PrognosticAnalysisEnum){ 361 362 UpdateInputsPrognostic( solution,analysis_type,sub_analysis_type); 363 } 364 else if (analysis_type==Prognostic2AnalysisEnum){ 365 366 UpdateInputsPrognostic2(solution,analysis_type,sub_analysis_type); 367 } 368 else if (analysis_type==BalancedthicknessAnalysisEnum){ 369 370 UpdateInputsBalancedthickness( solution,analysis_type,sub_analysis_type); 371 } 372 else if (analysis_type==Balancedthickness2AnalysisEnum){ 373 374 UpdateInputsBalancedthickness2( solution,analysis_type,sub_analysis_type); 375 } 376 else if (analysis_type==BalancedvelocitiesAnalysisEnum){ 377 378 UpdateInputsBalancedvelocities( solution,analysis_type,sub_analysis_type); 379 } 380 else{ 381 382 ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet"); 383 } 384 } 385 /*Object functions*/ 386 /*FUNCTION Penta::UpdateInputsDiagnosticHoriz {{{1*/ 387 void Penta::UpdateInputsDiagnosticHoriz(double* solution, int analysis_type, int sub_analysis_type){ 388 389 390 int i; 391 392 const int numvertices=6; 393 const int numdofpervertex=2; 394 const int numdof=numdofpervertex*numvertices; 395 396 int doflist[numdof]; 397 double values[numdof]; 398 double vx[numvertices]; 399 double vy[numvertices]; 400 401 int dummy; 402 403 /*Get dof list: */ 404 GetDofList(&doflist[0],&dummy); 405 406 /*Use the dof list to index into the solution vector: */ 407 for(i=0;i<numdof;i++){ 408 values[i]=solution[doflist[i]]; 409 } 410 411 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ 412 for(i=0;i<numvertices;i++){ 413 vx[i]=values[i*numdofpervertex+0]; 414 vy[i]=values[i*numdofpervertex+1]; 415 } 416 417 /*Add vx and vy as inputs to the tria element: */ 418 this->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); 419 this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); 420 } 421 422 /*}}}*/ 423 /*FUNCTION Penta::UpdateInputsSlopeCompute {{{1*/ 424 void Penta::UpdateInputsSlopeCompute(double* solution, int analysis_type, int sub_analysis_type){ 425 ISSMERROR(" not supported yet!"); 426 } 427 /*}}}*/ 428 /*FUNCTION Penta::UpdateInputsPrognostic {{{1*/ 429 void Penta::UpdateInputsPrognostic(double* solution, int analysis_type, int sub_analysis_type){ 430 ISSMERROR(" not supported yet!"); 431 } 432 /*}}}*/ 433 /*FUNCTION Penta::UpdateInputsPrognostic2 {{{1*/ 434 void Penta::UpdateInputsPrognostic2(double* solution, int analysis_type, int sub_analysis_type){ 435 ISSMERROR(" not supported yet!"); 436 } 437 /*}}}*/ 438 /*FUNCTION Penta::UpdateInputsBalancedthickness {{{1*/ 439 void Penta::UpdateInputsBalancedthickness(double* solution, int analysis_type, int sub_analysis_type){ 440 ISSMERROR(" not supported yet!"); 441 } 442 /*}}}*/ 443 /*FUNCTION Penta::UpdateInputsBalancedthickness2 {{{1*/ 444 void Penta::UpdateInputsBalancedthickness2(double* solution, int analysis_type, int sub_analysis_type){ 445 ISSMERROR(" not supported yet!"); 446 } 447 /*}}}*/ 448 /*FUNCTION Penta::UpdateInputsBalancedvelocities {{{1*/ 449 void Penta::UpdateInputsBalancedvelocities(double* solution, int analysis_type, int sub_analysis_type){ 277 450 ISSMERROR(" not supported yet!"); 278 451 } … … 295 468 double bed; 296 469 double basalforce[3]; 297 double vxvyvz_list[numgrids][3];298 double pressure_list[numgrids];299 470 double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ 300 471 double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ … … 304 475 int dofv[3]={0,1,2}; 305 476 int dofp[1]={3}; 306 ParameterInputs* inputs=NULL;307 477 double Jdet2d; 308 478 Tria* tria=NULL; … … 323 493 double surface=0; 324 494 double value=0; 495 496 /*flags: */ 497 bool onbed; 325 498 326 499 /*dynamic objects pointed to by hooks: */ … … 333 506 if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!"); 334 507 335 /*recover pointers: */336 inputs=(ParameterInputs*)vinputs;337 338 508 /*recover objects from hooks: */ 339 509 nodes=(Node**)hnodes.deliverp(); … … 342 512 numpar=(Numpar*)hnumpar.delivers(); 343 513 344 if(!this->properties.onbed){ 514 /*recover some inputs: */ 515 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 516 517 if(!onbed){ 345 518 //put zero 346 519 VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES); … … 351 524 rho_ice=matpar->GetRhoIce(); 352 525 gravity=matpar->GetG(); 353 354 /*recover extra inputs from users, at current convergence iteration: */355 inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofv,numgrids,(void**)nodes);356 inputs->Recover("velocity",&pressure_list[0] ,1,dofp,numgrids,(void**)nodes);357 526 358 527 /* Get node coordinates and dof list: */ … … 378 547 379 548 /*Compute strain rate viscosity and pressure: */ 380 GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);549 inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); 381 550 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 382 GetParameterValue(&pressure,&pressure_list[0],&gauss_coord[0]);551 inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum); 383 552 384 553 /*Compute Stress*/ … … 417 586 int i; 418 587 const int numgrids=6; 419 int doflist[numgrids];588 int doflist[numgrids]; 420 589 double pressure[numgrids]; 421 590 double rho_ice,g; 422 double xyz_list[numgrids][3]; 591 double surface[numgrids]; 592 double xyz_list[numgrids][3]; 593 double gauss[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 594 595 /*inputs: */ 596 bool onwater; 423 597 424 598 /*dynamic objects pointed to by hooks: */ … … 430 604 matpar=(Matpar*)hmatpar.delivers(); 431 605 606 /*retrieve inputs :*/ 607 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 608 432 609 /*If on water, skip: */ 433 if( this->properties.onwater)return;610 if(onwater)return; 434 611 435 612 /*Get node data: */ … … 441 618 /*Get dof list on which we will plug the pressure values: */ 442 619 GetDofList1(&doflist[0]); 620 621 /*recover value of surface at grids: */ 622 inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); 443 623 444 624 /*pressure is lithostatic: */ … … 446 626 g=matpar->GetG(); 447 627 for(i=0;i<numgrids;i++){ 448 pressure[i]=rho_ice*g*( this->properties.s[i]-xyz_list[i][2]);628 pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); 449 629 } 450 630 … … 467 647 Tria* tria=NULL; 468 648 649 /*flags: */ 650 bool onbed; 651 bool onwater; 652 bool collapse; 653 bool onsurface; 654 655 /*recover some inputs: */ 656 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 657 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 658 inputs->GetParameterValue(&collapse,CollapseEnum); 659 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 660 469 661 /*If on water, return 0: */ 470 if( this->properties.onwater)return 0;662 if(onwater)return 0; 471 663 472 664 /*Bail out if this element if: 473 665 * -> Non collapsed and not on the surface 474 666 * -> collapsed (2d model) and not on bed) */ 475 if ((! this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){667 if ((!collapse && !onsurface) || (collapse && !onbed)){ 476 668 return 0; 477 669 } 478 else if ( this->properties.collapse){670 else if (collapse){ 479 671 480 672 /*This element should be collapsed into a tria element at its base. Create this tria element, 481 673 * and compute CostFunction*/ 482 674 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 483 J=tria->CostFunction( inputs,analysis_type,sub_analysis_type);675 J=tria->CostFunction(analysis_type,sub_analysis_type); 484 676 delete tria; 485 677 return J; … … 488 680 489 681 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 490 J=tria->CostFunction( inputs,analysis_type,sub_analysis_type);682 J=tria->CostFunction(analysis_type,sub_analysis_type); 491 683 delete tria; 492 684 return J; … … 501 693 if (analysis_type==ControlAnalysisEnum){ 502 694 503 CreateKMatrixDiagnosticHoriz( Kgg, inputs,analysis_type,sub_analysis_type);695 CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type); 504 696 505 697 } … … 508 700 if (sub_analysis_type==HorizAnalysisEnum){ 509 701 510 CreateKMatrixDiagnosticHoriz( Kgg, inputs,analysis_type,sub_analysis_type);702 CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type); 511 703 } 512 704 else if (sub_analysis_type==VertAnalysisEnum){ 513 705 514 CreateKMatrixDiagnosticVert( Kgg, inputs,analysis_type,sub_analysis_type);706 CreateKMatrixDiagnosticVert( Kgg,analysis_type,sub_analysis_type); 515 707 } 516 708 else if (sub_analysis_type==StokesAnalysisEnum){ 517 709 518 CreateKMatrixDiagnosticStokes( Kgg, inputs,analysis_type,sub_analysis_type);710 CreateKMatrixDiagnosticStokes( Kgg,analysis_type,sub_analysis_type); 519 711 520 712 } … … 523 715 else if (analysis_type==SlopecomputeAnalysisEnum){ 524 716 525 CreateKMatrixSlopeCompute( Kgg, inputs,analysis_type,sub_analysis_type);717 CreateKMatrixSlopeCompute( Kgg,analysis_type,sub_analysis_type); 526 718 } 527 719 else if (analysis_type==PrognosticAnalysisEnum){ 528 720 529 CreateKMatrixPrognostic( Kgg, inputs,analysis_type,sub_analysis_type);721 CreateKMatrixPrognostic( Kgg,analysis_type,sub_analysis_type); 530 722 } 531 723 else if (analysis_type==BalancedthicknessAnalysisEnum){ … … 539 731 else if (analysis_type==ThermalAnalysisEnum){ 540 732 541 CreateKMatrixThermal( Kgg, inputs,analysis_type,sub_analysis_type);733 CreateKMatrixThermal( Kgg,analysis_type,sub_analysis_type); 542 734 } 543 735 else if (analysis_type==MeltingAnalysisEnum){ 544 736 545 CreateKMatrixMelting( Kgg, inputs,analysis_type,sub_analysis_type);737 CreateKMatrixMelting( Kgg,analysis_type,sub_analysis_type); 546 738 } 547 739 else{ … … 645 837 double Bprime[5][numdof]; 646 838 double L[2][numdof]; 647 double D[5][5]={ { 0,0,0,0,0 },{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}};// material matrix, simple scalar matrix.839 double D[5][5]={0.0}; // material matrix, simple scalar matrix. 648 840 double D_scalar; 649 double DL[2][2]={ { 0,0 },{0,0}}; //for basal drag841 double DL[2][2]={0.0}; //for basal drag 650 842 double DL_scalar; 651 843 652 844 /* local element matrices: */ 653 double Ke_gg[numdof][numdof]; //local element stiffness matrix 845 double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 846 654 847 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. 655 848 double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag … … 657 850 658 851 /*slope: */ 659 double slope[2]={0.0 ,0.0};852 double slope[2]={0.0}; 660 853 double slope_magnitude; 661 662 /*input parameters for structural analysis (diagnostic): */663 double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};664 double oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};665 int dofs[2]={0,1};666 double thickness;667 854 668 855 /*friction: */ … … 672 859 double MAXSLOPE=.06; // 6 % 673 860 double MOUNTAINKEXPONENT=10; 674 675 ParameterInputs* inputs=NULL;676 861 677 862 /*Collapsed formulation: */ … … 688 873 numpar=(Numpar*)hnumpar.delivers(); 689 874 875 /*inputs: */ 876 bool onwater; 877 bool collapse; 878 bool onbed; 879 bool shelf; 880 881 /*retrieve inputs :*/ 882 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 883 inputs->GetParameterValue(&collapse,CollapseEnum); 884 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 885 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 886 690 887 /*If on water, skip stiffness: */ 691 if(this->properties.onwater)return; 692 693 /*recover pointers: */ 694 inputs=(ParameterInputs*)vinputs; 888 if(onwater)return; 695 889 696 890 /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the … … 699 893 700 894 701 if (( this->properties.collapse==1) && (this->properties.onbed==0)){895 if ((collapse==1) && (onbed==0)){ 702 896 /*This element should be collapsed, but this element is not on the bedrock, therefore all its 703 897 * dofs have already been frozen! Do nothing: */ 704 898 return; 705 899 } 706 else if (( this->properties.collapse==1) && (this->properties.onbed==1)){900 else if ((collapse==1) && (onbed==1)){ 707 901 708 902 /*This element should be collapsed into a tria element at its base. Create this tria element, 709 903 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ 710 904 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 711 tria->CreateKMatrix(Kgg, inputs,analysis_type,sub_analysis_type);905 tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); 712 906 delete tria; 713 907 return; … … 716 910 717 911 /*Implement standard penta element: */ 718 719 /*recover extra inputs from users, at current convergence iteration: */720 inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);721 inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);722 912 723 913 /* Get node coordinates and dof list: */ … … 725 915 GetDofList(&doflist[0],&numberofdofspernode); 726 916 727 /* Set Ke_gg to 0: */728 for(i=0;i<numdof;i++){729 for(j=0;j<numdof;j++){730 Ke_gg[i][j]=0.0;731 }732 }733 917 734 918 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 760 944 761 945 /*Get strain rate from velocity: */ 762 GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_coord);763 GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_coord);946 inputs->GetStrainRate(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum); 947 inputs->GetStrainRate(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum); 764 948 765 949 /*Get viscosity: */ … … 774 958 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); 775 959 776 /*Build the D matrix: we plug the gaussian weight, the thickness, theviscosity, and the jacobian determinant960 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 777 961 onto this scalar matrix, so that we win some computational time: */ 778 962 … … 803 987 804 988 //Deal with 2d friction at the bedrock interface 805 if(( this->properties.onbed && !this->properties.shelf)){989 if((onbed && !shelf)){ 806 990 807 991 /*Build a tria element using the 3 grids of the base of the penta. Then use … … 810 994 811 995 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 812 tria->CreateKMatrixDiagnosticHorizFriction(Kgg, inputs,analysis_type,sub_analysis_type);996 tria->CreateKMatrixDiagnosticHorizFriction(Kgg,analysis_type,sub_analysis_type); 813 997 delete tria; 814 998 } … … 845 1029 int dofs[3]={0,1,2}; 846 1030 847 double K_terms[numdof][numdof] ;1031 double K_terms[numdof][numdof]={0.0}; 848 1032 849 1033 /*Material properties: */ … … 860 1044 double surface_normal[3]; 861 1045 double bed_normal[3]; 862 double vxvyvz_list[numgrids][3];863 1046 double thickness; 864 1047 … … 889 1072 double* area_gauss_weights = NULL; 890 1073 double* vert_gauss_weights = NULL; 1074 double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 891 1075 892 1076 /* specific gaussian point: */ … … 902 1086 double alpha2_gauss; 903 1087 904 ParameterInputs* inputs=NULL; 1088 double vx_list[numgrids]; 1089 double vy_list[numgrids]; 1090 double vz_list[numgrids]; 1091 double thickness_list[numgrids]; 1092 double bed_list[numgrids]; 1093 double dragcoefficient_list[numgrids]; 1094 double drag_p,drag_q; 905 1095 906 1096 /*dynamic objects pointed to by hooks: */ … … 910 1100 Numpar* numpar=NULL; 911 1101 1102 /*inputs: */ 1103 bool onwater; 1104 bool onbed; 1105 bool shelf; 1106 1107 /*retrieve inputs :*/ 1108 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1109 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1110 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 1111 912 1112 /*If on water, skip stiffness: */ 913 if( this->properties.onwater)return;1113 if(onwater)return; 914 1114 915 1115 /*recover objects from hooks: */ … … 919 1119 numpar=(Numpar*)hnumpar.delivers(); 920 1120 921 /*recover pointers: */922 inputs=(ParameterInputs*)vinputs;923 924 /* Set K_terms to 0: */925 for(i=0;i<numdof;i++){926 for(j=0;j<numdof;j++){927 K_terms[i][j]=0.0;928 }929 }930 1121 931 1122 /*recovre material parameters: */ … … 937 1128 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 938 1129 GetDofList(&doflist[0],&numberofdofspernode); 939 940 /*recover extra inputs from users, at current convergence iteration: */941 inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);942 1130 943 1131 /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 965 1153 gauss_coord[3]=*(vert_gauss_coord+igvert); 966 1154 967 /*Compute thickness: */968 969 1155 /*Compute strain rate: */ 970 GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);1156 inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); 971 1157 972 1158 /*Get viscosity: */ … … 1003 1189 } 1004 1190 1005 if(( this->properties.onbed==1) && (this->properties.shelf==0)){1191 if((onbed==1) && (shelf==0)){ 1006 1192 1007 1193 /*Build alpha2_list used by drag stiffness matrix*/ … … 1012 1198 strcpy(friction->element_type,"2d"); 1013 1199 1200 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],6,VxEnum); 1201 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],6,VyEnum); 1202 inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],6,VzEnum); 1203 inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],6,DragCoefficientEnum); 1204 inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],6,BedEnum); 1205 inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],6,ThicknessEnum); 1206 inputs->GetParameterValue(&drag_p,DragPEnum); 1207 inputs->GetParameterValue(&drag_q,DragQEnum); 1208 1014 1209 friction->gravity=matpar->GetG(); 1015 1210 friction->rho_ice=matpar->GetRhoIce(); 1016 1211 friction->rho_water=matpar->GetRhoWater(); 1017 friction->K=&this->properties.k[0]; 1018 friction->bed=&this->properties.b[0]; 1019 friction->thickness=&this->properties.h[0]; 1020 friction->velocities=&vxvyvz_list[0][0]; 1021 friction->p=this->properties.p; 1022 friction->q=this->properties.q; 1023 friction->analysis_type=analysis_type; 1212 friction->K=&dragcoefficient_list[0]; 1213 friction->bed=&bed_list[0]; 1214 friction->thickness=&thickness_list[0]; 1215 friction->vx=&vx_list[0]; 1216 friction->vy=&vy_list[0]; 1217 friction->vz=&vz_list[0]; 1218 friction->p=drag_p; 1219 friction->q=drag_q; 1024 1220 1025 1221 /*Compute alpha2_list: */ … … 1058 1254 1059 1255 /*Compute strain rate: */ 1060 GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);1256 inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); 1061 1257 1062 1258 /*Get viscosity at last iteration: */ … … 1098 1294 } 1099 1295 } 1100 } //if ( ( this->properties.onbed==1) && (shelf==0))1296 } //if ( (onbed==1) && (shelf==0)) 1101 1297 1102 1298 /*Reduce the matrix */ … … 1155 1351 1156 1352 /* matrices: */ 1157 double Ke_gg[numdof][numdof] ;1353 double Ke_gg[numdof][numdof]={0.0}; 1158 1354 double Ke_gg_gaussian[numdof][numdof]; 1159 1355 double Jdet; … … 1168 1364 nodes=(Node**)hnodes.deliverp(); 1169 1365 1170 ParameterInputs* inputs=NULL;1171 1172 1366 /*Collapsed formulation: */ 1173 1367 Tria* tria=NULL; 1174 1368 1369 /*inputs: */ 1370 bool onwater; 1371 bool onsurface; 1372 1373 /*retrieve inputs :*/ 1374 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1375 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 1376 1175 1377 /*If on water, skip stiffness: */ 1176 if(this->properties.onwater)return; 1177 1178 /*recover pointers: */ 1179 inputs=(ParameterInputs*)vinputs; 1180 1378 if(onwater)return; 1181 1379 1182 1380 /*If this element is on the surface, we have a dynamic boundary condition that applies, as a stiffness 1183 1381 * matrix: */ 1184 if( this->properties.onsurface){1382 if(onsurface){ 1185 1383 tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface 1186 tria->CreateKMatrixDiagnosticSurfaceVert(Kgg, inputs,analysis_type,sub_analysis_type);1384 tria->CreateKMatrixDiagnosticSurfaceVert(Kgg, analysis_type,sub_analysis_type); 1187 1385 delete tria; 1188 1386 } … … 1194 1392 GetDofList(&doflist[0],&numberofdofspernode); 1195 1393 1196 /* Set Ke_gg to 0: */1197 for(i=0;i<numdof;i++){1198 for(j=0;j<numdof;j++){1199 Ke_gg[i][j]=0.0;1200 }1201 }1202 1394 1203 1395 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 1265 1457 Tria* tria=NULL; 1266 1458 1459 /*inputs: */ 1460 bool onwater; 1461 bool onbed; 1462 1463 /*retrieve inputs :*/ 1464 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1465 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1466 1267 1467 /*If on water, skip: */ 1268 if( this->properties.onwater)return;1269 1270 if (! this->properties.onbed){1468 if(onwater)return; 1469 1470 if (!onbed){ 1271 1471 return; 1272 1472 } … … 1274 1474 1275 1475 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 1276 tria->CreateKMatrixMelting(Kgg, inputs,analysis_type,sub_analysis_type);1476 tria->CreateKMatrixMelting(Kgg, analysis_type,sub_analysis_type); 1277 1477 delete tria; 1278 1478 return; … … 1287 1487 Tria* tria=NULL; 1288 1488 1489 /*inputs: */ 1490 bool onwater; 1491 bool onbed; 1492 1493 /*retrieve inputs :*/ 1494 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1495 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1496 1289 1497 /*If on water, skip: */ 1290 if( this->properties.onwater)return;1498 if(onwater)return; 1291 1499 1292 1500 /*Is this element on the bed? :*/ 1293 if(! this->properties.onbed)return;1501 if(!onbed)return; 1294 1502 1295 1503 /*Spawn Tria element from the base of the Penta: */ 1296 1504 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 1297 tria->CreateKMatrix(Kgg, inputs,analysis_type,sub_analysis_type);1505 tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); 1298 1506 delete tria; 1299 1507 return; … … 1308 1516 Tria* tria=NULL; 1309 1517 1518 /*inputs: */ 1519 bool onwater; 1520 bool onbed; 1521 1522 /*retrieve inputs :*/ 1523 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1524 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1525 1526 1310 1527 /*If on water, skip: */ 1311 if( this->properties.onwater)return;1528 if(onwater)return; 1312 1529 1313 1530 /*Is this element on the bed? :*/ 1314 if(! this->properties.onbed)return;1531 if(!onbed)return; 1315 1532 1316 1533 /*Spawn Tria element from the base of the Penta: */ 1317 1534 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 1318 tria->CreateKMatrix(Kgg, inputs,analysis_type,sub_analysis_type);1535 tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); 1319 1536 delete tria; 1320 1537 return; … … 1356 1573 double K[2][2]={0.0}; 1357 1574 1358 double vxvyvz_list[numgrids][3]; 1359 double vx_list[numgrids]; 1360 int vx_list_exists; 1361 double u; 1362 double vy_list[numgrids]; 1363 int vy_list_exists; 1364 double v; 1365 double vz_list[numgrids]; 1366 int vz_list_exists; 1367 double w; 1575 double u,v,w; 1368 1576 1369 1577 /*matrices: */ … … 1404 1612 /*Collapsed formulation: */ 1405 1613 Tria* tria=NULL; 1406 ParameterInputs* inputs=NULL; 1614 1615 1616 /*inputs: */ 1617 bool onwater; 1618 bool onbed; 1619 bool shelf; 1620 1621 /*retrieve inputs :*/ 1622 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1623 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1624 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 1407 1625 1408 1626 /*If on water, skip: */ 1409 if( this->properties.onwater)return;1627 if(onwater)return; 1410 1628 1411 1629 /*recover objects from hooks: */ … … 1413 1631 matpar=(Matpar*)hmatpar.delivers(); 1414 1632 numpar=(Numpar*)hnumpar.delivers(); 1415 1416 /*recover pointers: */1417 inputs=(ParameterInputs*)vinputs;1418 1633 1419 1634 /* Get node coordinates and dof list: */ … … 1427 1642 heatcapacity=matpar->GetHeatCapacity(); 1428 1643 thermalconductivity=matpar->GetThermalConductivity(); 1429 1430 /*recover extra inputs from users, dt and velocity: */ 1431 found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); 1432 if(!found)ISSMERROR(" could not find velocity in inputs!"); 1433 found=inputs->Recover("dt",&dt); 1434 if(!found)ISSMERROR(" could not find dt in inputs!"); 1435 1436 for(i=0;i<numgrids;i++){ 1437 vx_list[i]=vxvyvz_list[i][0]; 1438 vy_list[i]=vxvyvz_list[i][1]; 1439 vz_list[i]=vxvyvz_list[i][2]; 1440 } 1441 1644 dt=numpar->dt; 1442 1645 1443 1646 /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 1495 1698 1496 1699 //Build the D matrix 1497 GetParameterValue(&u, &vx_list[0],gauss_coord);1498 GetParameterValue(&v, &vy_list[0],gauss_coord);1499 GetParameterValue(&w, &vz_list[0],gauss_coord);1700 inputs->GetParameterValue(&u, gauss_coord,VxEnum); 1701 inputs->GetParameterValue(&v, gauss_coord,VyEnum); 1702 inputs->GetParameterValue(&w, gauss_coord,VzEnum); 1500 1703 1501 1704 D_scalar=gauss_weight*Jdet; … … 1577 1780 1578 1781 //Ice/ocean heat exchange flux on ice shelf base 1579 if( this->properties.onbed && this->properties.shelf){1782 if(onbed && shelf){ 1580 1783 1581 1784 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 1582 tria->CreateKMatrixThermal(Kgg, inputs,analysis_type,sub_analysis_type);1785 tria->CreateKMatrixThermal(Kgg, analysis_type,sub_analysis_type); 1583 1786 delete tria; 1584 1787 } … … 1591 1794 if (analysis_type==ControlAnalysisEnum){ 1592 1795 1593 CreatePVectorDiagnosticHoriz( pg, inputs,analysis_type,sub_analysis_type);1796 CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type); 1594 1797 } 1595 1798 else if (analysis_type==DiagnosticAnalysisEnum){ … … 1597 1800 if (sub_analysis_type==HorizAnalysisEnum){ 1598 1801 1599 CreatePVectorDiagnosticHoriz( pg, inputs,analysis_type,sub_analysis_type);1802 CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type); 1600 1803 } 1601 1804 else if (sub_analysis_type==VertAnalysisEnum){ 1602 1805 1603 CreatePVectorDiagnosticVert( pg, inputs,analysis_type,sub_analysis_type);1806 CreatePVectorDiagnosticVert( pg,analysis_type,sub_analysis_type); 1604 1807 } 1605 1808 else if (sub_analysis_type==StokesAnalysisEnum){ 1606 1809 1607 CreatePVectorDiagnosticStokes( pg, inputs,analysis_type,sub_analysis_type);1810 CreatePVectorDiagnosticStokes( pg,analysis_type,sub_analysis_type); 1608 1811 } 1609 1812 else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); … … 1611 1814 else if (analysis_type==SlopecomputeAnalysisEnum){ 1612 1815 1613 CreatePVectorSlopeCompute( pg, inputs,analysis_type,sub_analysis_type);1816 CreatePVectorSlopeCompute( pg,analysis_type,sub_analysis_type); 1614 1817 } 1615 1818 else if (analysis_type==PrognosticAnalysisEnum){ 1616 1819 1820 CreatePVectorPrognostic( pg,analysis_type,sub_analysis_type); 1821 } 1822 else if (analysis_type==BalancedthicknessAnalysisEnum){ 1823 1617 1824 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); 1618 1825 } 1619 else if (analysis_type==Balanced thicknessAnalysisEnum){1826 else if (analysis_type==BalancedvelocitiesAnalysisEnum){ 1620 1827 1621 1828 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); 1622 1829 } 1623 else if (analysis_type==BalancedvelocitiesAnalysisEnum){1624 1625 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);1626 }1627 1830 else if (analysis_type==ThermalAnalysisEnum){ 1628 1831 1629 CreatePVectorThermal( pg, inputs,analysis_type,sub_analysis_type);1832 CreatePVectorThermal( pg,analysis_type,sub_analysis_type); 1630 1833 } 1631 1834 else if (analysis_type==MeltingAnalysisEnum){ 1632 1835 1633 CreatePVectorMelting( pg, inputs,analysis_type,sub_analysis_type);1836 CreatePVectorMelting( pg,analysis_type,sub_analysis_type); 1634 1837 } 1635 1838 else{ … … 1720 1923 1721 1924 /*element vector at the gaussian points: */ 1722 double pe_g[numdof] ;1925 double pe_g[numdof]={0.0}; 1723 1926 double pe_g_gaussian[numdof]; 1724 1725 ParameterInputs* inputs=NULL;1726 1927 1727 1928 /*dynamic objects pointed to by hooks: */ … … 1732 1933 Tria* tria=NULL; 1733 1934 1935 /*inputs: */ 1936 bool onwater; 1937 bool collapse; 1938 bool onbed; 1939 1940 /*retrieve inputs :*/ 1941 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1942 inputs->GetParameterValue(&collapse,CollapseEnum); 1943 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1944 1734 1945 /*If on water, skip load: */ 1735 if(this->properties.onwater)return; 1736 1737 /*recover pointers: */ 1738 inputs=(ParameterInputs*)vinputs; 1946 if(onwater)return; 1739 1947 1740 1948 /*recover objects from hooks: */ … … 1747 1955 the load vector. */ 1748 1956 1749 if (( this->properties.collapse==1) && (this->properties.onbed==0)){1957 if ((collapse==1) && (onbed==0)){ 1750 1958 /*This element should be collapsed, but this element is not on the bedrock, therefore all its 1751 1959 * dofs have already been frozen! Do nothing: */ 1752 1960 return; 1753 1961 } 1754 else if (( this->properties.collapse==1) && (this->properties.onbed==1)){1962 else if ((collapse==1) && (onbed==1)){ 1755 1963 1756 1964 /*This element should be collapsed into a tria element at its base. Create this tria element, 1757 1965 *and use its CreatePVector functionality to return an elementary load vector: */ 1758 1966 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 1759 tria->CreatePVector(pg, inputs,analysis_type,sub_analysis_type);1967 tria->CreatePVector(pg, analysis_type,sub_analysis_type); 1760 1968 delete tria; 1761 1969 return; … … 1768 1976 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 1769 1977 GetDofList(&doflist[0],&numberofdofspernode); 1770 1771 /* Set pe_g to 0: */1772 for(i=0;i<numdof;i++) pe_g[i]=0.0;1773 1978 1774 1979 /*Get gaussian points and weights :*/ … … 1793 1998 1794 1999 /*Compute thickness at gaussian point: */ 1795 GetParameterValue(&thickness, &this->properties.h[0],gauss_coord);2000 inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum); 1796 2001 1797 2002 /*Compute slope at gaussian point: */ 1798 GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_coord);2003 inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum); 1799 2004 1800 2005 /* Get Jacobian determinant: */ … … 1820 2025 } //for (ig1=0; ig1<num_area_gauss; ig1++) 1821 2026 1822 } //else if ((collapse==1) && ( this->properties.onbed==1))2027 } //else if ((collapse==1) && (onbed==1)) 1823 2028 1824 2029 /*Add pe_g to global vector pg: */ … … 1857 2062 double bed_normal[3]; 1858 2063 double bed; 1859 double vxvyvz_list[numgrids][3];1860 2064 1861 2065 /* gaussian points: */ … … 1898 2102 double D_scalar; 1899 2103 double tBD[27][8]; 1900 double P_terms[numdof]; 1901 1902 ParameterInputs* inputs=NULL; 2104 double P_terms[numdof]={0.0}; 2105 1903 2106 Tria* tria=NULL; 1904 2107 … … 1909 2112 Numpar* numpar=NULL; 1910 2113 2114 /*inputs: */ 2115 bool onwater; 2116 bool onbed; 2117 bool shelf; 2118 2119 /*retrieve inputs :*/ 2120 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2121 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2122 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 2123 1911 2124 /*If on water, skip load: */ 1912 if(this->properties.onwater)return; 1913 1914 /*recover pointers: */ 1915 inputs=(ParameterInputs*)vinputs; 2125 if(onwater)return; 1916 2126 1917 2127 /*recover objects from hooks: */ … … 1921 2131 numpar=(Numpar*)hnumpar.delivers(); 1922 2132 1923 /* Set P_terms to 0: */1924 for(i=0;i<numdof;i++){1925 P_terms[i]=0;1926 }1927 2133 1928 2134 /*recovre material parameters: */ … … 1930 2136 rho_ice=matpar->GetRhoIce(); 1931 2137 gravity=matpar->GetG(); 1932 1933 /*recover extra inputs from users, at current convergence iteration: */1934 inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);1935 2138 1936 2139 /* Get node coordinates and dof list: */ … … 1960 2163 1961 2164 /*Compute strain rate and viscosity: */ 1962 GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);2165 inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); 1963 2166 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 1964 2167 … … 2014 2217 2015 2218 /*Deal with 2d friction at the bedrock interface: */ 2016 if ( ( this->properties.onbed==1) && (this->properties.shelf==1)){2219 if ( (onbed==1) && (shelf==1)){ 2017 2220 2018 2221 for(i=0;i<numgrids2d;i++){ … … 2041 2244 2042 2245 /* Get bed at gaussian point */ 2043 GetParameterValue(&bed,&this->properties.b[0],gauss_coord);2246 inputs->GetParameterValue(&bed, gauss_coord,BedEnum); 2044 2247 2045 2248 /*Get L matrix : */ … … 2062 2265 } 2063 2266 } 2064 } //if ( ( this->properties.onbed==1) && (this->properties.shelf==1))2267 } //if ( (onbed==1) && (shelf==1)) 2065 2268 2066 2269 /*Reduce the matrix */ … … 2118 2321 2119 2322 /*element vector at the gaussian points: */ 2120 double pe_g[numdof] ;2323 double pe_g[numdof]={0.0}; 2121 2324 double pe_g_gaussian[numdof]; 2122 2325 double l1l6[6]; … … 2126 2329 2127 2330 /*input parameters for structural analysis (diagnostic): */ 2128 double vx_list[numgrids]={0,0,0,0,0,0};2129 double vy_list[numgrids]={0,0,0,0,0,0};2130 2331 double du[3]; 2131 2332 double dv[3]; … … 2137 2338 Node** nodes=NULL; 2138 2339 2139 ParameterInputs* inputs=NULL; 2340 /*inputs: */ 2341 bool onwater; 2342 bool onbed; 2140 2343 2141 2344 /*recover objects from hooks: */ 2142 2345 nodes=(Node**)hnodes.deliverp(); 2143 2346 2144 /*recover pointers: */ 2145 inputs=(ParameterInputs*)vinputs; 2146 2347 /*retrieve inputs :*/ 2348 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2349 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2350 2351 2147 2352 /*If on water, skip: */ 2148 if( this->properties.onwater)return;2353 if(onwater)return; 2149 2354 2150 2355 /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the 2151 2356 *diagnostic base vertical stifness: */ 2152 if( this->properties.onbed){2357 if(onbed){ 2153 2358 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock 2154 tria->CreatePVectorDiagnosticBaseVert(pg, inputs,analysis_type,sub_analysis_type);2359 tria->CreatePVectorDiagnosticBaseVert(pg, analysis_type,sub_analysis_type); 2155 2360 delete tria; 2156 2361 } … … 2161 2366 GetDofList(&doflist[0],&numberofdofspernode); 2162 2367 2163 /* Set pe_g to 0: */2164 for(i=0;i<numdof;i++) pe_g[i]=0.0;2165 2166 /* recover input parameters: */2167 if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))ISSMERROR(" cannot compute vertical velocity without horizontal velocity!");2168 inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);2169 2170 2368 /*Get gaussian points and weights :*/ 2171 2369 order_area_gauss=2; … … 2189 2387 2190 2388 /*Get velocity derivative, with respect to x and y: */ 2191 GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_coord); 2192 GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_coord); 2389 2390 inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum); 2391 inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VyEnum); 2193 2392 dudx=du[0]; 2194 2393 dvdy=dv[1]; … … 2235 2434 Tria* tria=NULL; 2236 2435 2436 /*inputs: */ 2437 bool onwater; 2438 bool onbed; 2439 2440 /*retrieve inputs :*/ 2441 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2442 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2443 2237 2444 /*If on water, skip: */ 2238 if( this->properties.onwater)return;2445 if(onwater)return; 2239 2446 2240 2447 /*Is this element on the bed? :*/ 2241 if(! this->properties.onbed)return;2448 if(!onbed)return; 2242 2449 2243 2450 /*Spawn Tria element from the base of the Penta: */ 2244 2451 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2245 tria->CreatePVector(pg, inputs,analysis_type,sub_analysis_type);2452 tria->CreatePVector(pg, analysis_type,sub_analysis_type); 2246 2453 delete tria; 2247 2454 return; … … 2255 2462 Tria* tria=NULL; 2256 2463 2464 /*inputs: */ 2465 bool onwater; 2466 bool onbed; 2467 2468 /*retrieve inputs :*/ 2469 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2470 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2471 2257 2472 /*If on water, skip: */ 2258 if( this->properties.onwater)return;2473 if(onwater)return; 2259 2474 2260 2475 /*Is this element on the bed? :*/ 2261 if(! this->properties.onbed)return;2476 if(!onbed)return; 2262 2477 2263 2478 /*Spawn Tria element from the base of the Penta: */ 2264 2479 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2265 tria->CreatePVector(pg, inputs,analysis_type,sub_analysis_type);2480 tria->CreatePVector(pg, analysis_type,sub_analysis_type); 2266 2481 delete tria; 2267 2482 return; … … 2299 2514 2300 2515 double dt; 2301 double vx_list[numgrids];2302 double vy_list[numgrids];2303 double vz_list[numgrids];2304 double vxvyvz_list[numgrids][3];2305 2516 double temperature_list[numgrids]; 2306 2517 double temperature; … … 2340 2551 /*Collapsed formulation: */ 2341 2552 Tria* tria=NULL; 2342 ParameterInputs* inputs=NULL;2343 2553 2344 2554 /*dynamic objects pointed to by hooks: */ … … 2346 2556 Matpar* matpar=NULL; 2347 2557 Matice* matice=NULL; 2558 Numpar* numpar=NULL; 2559 2560 /*inputs: */ 2561 bool onwater; 2562 bool onbed; 2563 bool shelf; 2564 2565 /*retrieve inputs :*/ 2566 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2567 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2568 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 2348 2569 2349 2570 /*If on water, skip: */ 2350 if(this->properties.onwater)return; 2351 2352 /*recover pointers: */ 2353 inputs=(ParameterInputs*)vinputs; 2571 if(onwater)return; 2354 2572 2355 2573 /*recover objects from hooks: */ … … 2357 2575 matpar=(Matpar*)hmatpar.delivers(); 2358 2576 matice=(Matice*)hmatice.delivers(); 2577 numpar=(Numpar*)hnumpar.delivers(); 2359 2578 2360 2579 /* Get node coordinates and dof list: */ … … 2369 2588 beta=matpar->GetBeta(); 2370 2589 meltingpoint=matpar->GetMeltingPoint(); 2371 2372 /*recover extra inputs from users, dt and velocity: */ 2373 found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); 2374 if(!found)ISSMERROR(" could not find velocity in inputs!"); 2375 found=inputs->Recover("dt",&dt); 2376 if(!found)ISSMERROR(" could not find dt in inputs!"); 2377 2378 if(dt){ 2379 found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes); 2380 if(!found)ISSMERROR(" could not find temperature in inputs!"); 2381 } 2382 2383 for(i=0;i<numgrids;i++){ 2384 vx_list[i]=vxvyvz_list[i][0]; 2385 vy_list[i]=vxvyvz_list[i][1]; 2386 vz_list[i]=vxvyvz_list[i][2]; 2387 } 2590 dt=numpar->dt; 2388 2591 2389 2592 /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore … … 2409 2612 2410 2613 /*Compute strain rate and viscosity: */ 2411 GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);2614 inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); 2412 2615 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 2413 2616 … … 2433 2636 /* Build transient now */ 2434 2637 if(dt){ 2435 GetParameterValue(&temperature, &temperature_list[0],gauss_coord);2638 inputs->GetParameterValue(&temperature, gauss_coord,TemperatureEnum); 2436 2639 scalar_transient=temperature*Jdet*gauss_weight; 2437 2640 for(i=0;i<numgrids;i++){ … … 2446 2649 2447 2650 /* Ice/ocean heat exchange flux on ice shelf base */ 2448 if( this->properties.onbed && this->properties.shelf){2651 if(onbed && shelf){ 2449 2652 2450 2653 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2451 tria->CreatePVectorThermalShelf(pg, inputs,analysis_type,sub_analysis_type);2654 tria->CreatePVectorThermalShelf(pg, analysis_type,sub_analysis_type); 2452 2655 delete tria; 2453 2656 } 2454 2657 2455 2658 /* Geothermal flux on ice sheet base and basal friction */ 2456 if( this->properties.onbed && !this->properties.shelf){2659 if(onbed && !shelf){ 2457 2660 2458 2661 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2459 tria->CreatePVectorThermalSheet(pg, inputs,analysis_type,sub_analysis_type);2662 tria->CreatePVectorThermalSheet(pg, analysis_type,sub_analysis_type); 2460 2663 delete tria; 2461 2664 } … … 2477 2680 Tria* tria=NULL; 2478 2681 2682 /*inputs: */ 2683 bool onwater; 2684 bool collapse; 2685 bool onsurface; 2686 bool onbed; 2687 2688 /*retrieve inputs :*/ 2689 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2690 inputs->GetParameterValue(&collapse,CollapseEnum); 2691 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2692 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 2693 2479 2694 /*If on water, skip: */ 2480 if( this->properties.onwater)return;2695 if(onwater)return; 2481 2696 2482 2697 /*Bail out if this element if: 2483 2698 * -> Non collapsed and not on the surface 2484 2699 * -> collapsed (2d model) and not on bed) */ 2485 if ((! this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){2700 if ((!collapse && !onsurface) || (collapse && !onbed)){ 2486 2701 return; 2487 2702 } 2488 else if ( this->properties.collapse){2703 else if (collapse){ 2489 2704 2490 2705 /*This element should be collapsed into a tria element at its base. Create this tria element, 2491 2706 * and compute Du*/ 2492 2707 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 2493 tria->Du(du_g, inputs,analysis_type,sub_analysis_type);2708 tria->Du(du_g,analysis_type,sub_analysis_type); 2494 2709 delete tria; 2495 2710 return; … … 2498 2713 2499 2714 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 2500 tria->Du(du_g, inputs,analysis_type,sub_analysis_type);2715 tria->Du(du_g,analysis_type,sub_analysis_type); 2501 2716 delete tria; 2502 2717 return; 2503 2718 } 2504 }2505 /*}}}*/2506 /*FUNCTION Enum {{{1*/2507 int Penta::Enum(void){2508 2509 return PentaEnum;2510 2511 2719 } 2512 2720 /*}}}*/ … … 2527 2735 nodes=(Node**)hnodes.deliverp(); 2528 2736 2737 /*inputs: */ 2738 bool collapse; 2739 bool onbed; 2740 2741 /*retrieve inputs :*/ 2742 inputs->GetParameterValue(&collapse,CollapseEnum); 2743 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2744 2529 2745 /*Figure out if we should extrude for this element: */ 2530 2746 if (iscollapsed){ 2531 2747 /*From higher level, we are told to extrude only elements that have the collapse flag on: */ 2532 if ( this->properties.collapse)extrude=1;2748 if (collapse)extrude=1; 2533 2749 else extrude=0; 2534 2750 } … … 2540 2756 /*Now, extrusion starts from the bed on, so double check this element is on 2541 2757 * the bedrock: */ 2542 if( this->properties.onbed==0)extrude=0;2758 if(onbed==0)extrude=0; 2543 2759 2544 2760 /*Go on and extrude field: */ … … 2822 3038 /*}}}*/ 2823 3039 /*FUNCTION GetBedList {{{1*/ 2824 void Penta::GetBedList(double* bed_list){ 2825 2826 int i; 2827 for(i=0;i<6;i++)bed_list[i]=this->properties.b[i]; 3040 void Penta::GetBedList(double* bedlist){ 3041 3042 const int numgrids=6; 3043 double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 3044 3045 inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum); 2828 3046 2829 3047 } … … 3682 3900 /*FUNCTION GetOnBed {{{1*/ 3683 3901 bool Penta::GetOnBed(){ 3684 return this->properties.onbed; 3902 3903 bool onbed; 3904 3905 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 3906 3907 return onbed; 3685 3908 } 3686 3909 /*}}}*/ … … 3757 3980 /*}}}*/ 3758 3981 /*FUNCTION GetShelf {{{1*/ 3759 int Penta::GetShelf(){ 3760 return this->properties.shelf; 3982 bool Penta::GetShelf(){ 3983 bool onshelf; 3984 3985 /*retrieve inputs :*/ 3986 inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum); 3987 3988 return onshelf; 3761 3989 } 3762 3990 /*}}}*/ … … 3821 4049 void Penta::GetThicknessList(double* thickness_list){ 3822 4050 3823 int i; 3824 for(i=0;i<6;i++)thickness_list[i]=this->properties.h[i]; 4051 const int numgrids=6; 4052 double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; 4053 inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum); 4054 3825 4055 } 3826 4056 /*}}}*/ … … 3828 4058 void Penta::Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type){ 3829 4059 4060 /*inputs: */ 4061 bool onwater; 4062 4063 /*retrieve inputs :*/ 4064 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 4065 3830 4066 /*If on water, skip grad (=0): */ 3831 if( this->properties.onwater)return;4067 if(onwater)return; 3832 4068 3833 4069 if (strcmp(control_type,"drag")==0){ 3834 GradjDrag( grad_g, inputs,analysis_type,sub_analysis_type);4070 GradjDrag( grad_g,analysis_type,sub_analysis_type); 3835 4071 } 3836 4072 else if (strcmp(control_type,"B")==0){ 3837 GradjB( grad_g, inputs,analysis_type,sub_analysis_type);4073 GradjB( grad_g, analysis_type,sub_analysis_type); 3838 4074 } 3839 4075 else ISSMERROR("%s%s","control type not supported yet: ",control_type); … … 3845 4081 Tria* tria=NULL; 3846 4082 4083 /*inputs: */ 4084 bool onwater; 4085 bool onbed; 4086 bool shelf; 4087 4088 /*retrieve inputs :*/ 4089 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 4090 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 4091 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 4092 3847 4093 /*If on water, skip: */ 3848 if( this->properties.onwater)return;4094 if(onwater)return; 3849 4095 3850 4096 /*If on shelf, skip: */ 3851 if( this->properties.shelf)return;4097 if(shelf)return; 3852 4098 3853 4099 /*Bail out if this element does not touch the bedrock: */ 3854 if (! this->properties.onbed) return;4100 if (!onbed) return; 3855 4101 3856 4102 if (sub_analysis_type==HorizAnalysisEnum){ … … 3858 4104 /*MacAyeal or Pattyn*/ 3859 4105 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 3860 tria->GradjDrag( grad_g, inputs,analysis_type,sub_analysis_type);4106 tria->GradjDrag( grad_g,analysis_type,sub_analysis_type); 3861 4107 delete tria; 3862 4108 return; … … 3866 4112 /*Stokes*/ 3867 4113 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 3868 tria->GradjDragStokes( grad_g, inputs,analysis_type,sub_analysis_type);4114 tria->GradjDragStokes( grad_g,analysis_type,sub_analysis_type); 3869 4115 delete tria; 3870 4116 return; … … 3878 4124 Tria* tria=NULL; 3879 4125 4126 /*inputs: */ 4127 bool onwater; 4128 bool collapse; 4129 bool onbed; 4130 4131 /*retrieve inputs :*/ 4132 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 4133 inputs->GetParameterValue(&collapse,CollapseEnum); 4134 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 4135 3880 4136 /*If on water, skip: */ 3881 if( this->properties.onwater)return;3882 3883 if ( this->properties.collapse){4137 if(onwater)return; 4138 4139 if (collapse){ 3884 4140 /*Bail out element if collapsed (2d) and not on bed*/ 3885 if (! this->properties.onbed) return;4141 if (!onbed) return; 3886 4142 3887 4143 /*This element should be collapsed into a tria element at its base. Create this tria element, 3888 4144 * and compute gardj*/ 3889 4145 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 3890 tria->GradjB(grad_g, inputs,analysis_type,sub_analysis_type);4146 tria->GradjB(grad_g,analysis_type,sub_analysis_type); 3891 4147 delete tria; 3892 4148 return; … … 3895 4151 /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ 3896 4152 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 3897 tria->GradjB(grad_g, inputs,analysis_type,sub_analysis_type);4153 tria->GradjB(grad_g,analysis_type,sub_analysis_type); 3898 4154 delete tria; 3899 4155 return; … … 3912 4168 Tria* tria=NULL; 3913 4169 4170 /*inputs: */ 4171 bool onwater; 4172 bool collapse; 4173 bool onsurface; 4174 bool onbed; 4175 4176 /*retrieve inputs :*/ 4177 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 4178 inputs->GetParameterValue(&collapse,CollapseEnum); 4179 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 4180 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 4181 3914 4182 /*If on water, return 0: */ 3915 if( this->properties.onwater)return 0;4183 if(onwater)return 0; 3916 4184 3917 4185 /*Bail out if this element if: 3918 4186 * -> Non collapsed and not on the surface 3919 4187 * -> collapsed (2d model) and not on bed) */ 3920 if ((! this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){4188 if ((!collapse && !onsurface) || (collapse && !onbed)){ 3921 4189 return 0; 3922 4190 } 3923 else if ( this->properties.collapse){4191 else if (collapse){ 3924 4192 3925 4193 /*This element should be collapsed into a tria element at its base. Create this tria element, 3926 4194 * and compute Misfit*/ 3927 4195 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 3928 J=tria->Misfit( inputs,analysis_type,sub_analysis_type);4196 J=tria->Misfit(analysis_type,sub_analysis_type); 3929 4197 delete tria; 3930 4198 return J; … … 3933 4201 3934 4202 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 3935 J=tria->Misfit( inputs,analysis_type,sub_analysis_type);4203 J=tria->Misfit(analysis_type,sub_analysis_type); 3936 4204 delete tria; 3937 4205 return J; … … 4045 4313 Tria* tria=NULL; 4046 4314 4315 /*inputs: */ 4316 bool onwater; 4317 bool collapse; 4318 bool onsurface; 4319 bool onbed; 4320 4321 /*retrieve inputs :*/ 4322 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 4323 inputs->GetParameterValue(&collapse,CollapseEnum); 4324 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 4325 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 4326 4047 4327 /*If on water, return 0: */ 4048 if( this->properties.onwater)return 0;4328 if(onwater)return 0; 4049 4329 4050 4330 /*Bail out if this element if: 4051 4331 * -> Non collapsed and not on the surface 4052 4332 * -> collapsed (2d model) and not on bed) */ 4053 if ((! this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){4333 if ((!collapse && !onsurface) || (collapse && !onbed)){ 4054 4334 return 0; 4055 4335 } 4056 else if ( this->properties.collapse){4336 else if (collapse){ 4057 4337 4058 4338 /*This element should be collapsed into a tria element at its base. Create this tria element, 4059 4339 * and compute SurfaceArea*/ 4060 4340 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 4061 S=tria->SurfaceArea( inputs,analysis_type,sub_analysis_type);4341 S=tria->SurfaceArea(analysis_type,sub_analysis_type); 4062 4342 delete tria; 4063 4343 return S; … … 4066 4346 4067 4347 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 4068 S=tria->SurfaceArea( inputs,analysis_type,sub_analysis_type);4348 S=tria->SurfaceArea(analysis_type,sub_analysis_type); 4069 4349 delete tria; 4070 4350 return S; -
issm/trunk/src/c/objects/Penta.h
r3613 r3620 16 16 #include "./Element.h" 17 17 #include "./Matpar.h" 18 #include "./Numpar.h"19 18 #include "./Matice.h" 20 19 #include "./Tria.h" … … 45 44 /*}}}*/ 46 45 /*FUNCTION object management {{{1*/ 47 void Configure( void* loads,void* nodes,void* materials,void* parameters);46 void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); 48 47 Object* copy(); 49 48 void DeepEcho(); … … 57 56 int MyRank(); 58 57 void* SpawnTria(int g0, int g1, int g2); 59 void UpdateFromDakota(void* inputs);60 void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type);61 58 void SetClone(int* minranks); 62 59 … … 107 104 void CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type); 108 105 void CreatePVectorPrognostic( Vec pg, int analysis_type,int sub_analysis_type); 109 void CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);110 void CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);111 void CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);112 void CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);113 106 114 107 void CreateKMatrixDiagnosticStokes( Mat Kgg, int analysis_type,int sub_analysis_type); … … 136 129 void GetPhi(double* phi, double* epsilon, double viscosity); 137 130 double MassFlux(double* segment,double* ug); 131 132 /*updates: */ 133 void UpdateFromDakota(void* inputs); 134 void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); 135 void UpdateInputsDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type); 136 void UpdateInputsSlopeCompute( double* solution,int analysis_type,int sub_analysis_type); 137 void UpdateInputsPrognostic( double* solution,int analysis_type,int sub_analysis_type); 138 void UpdateInputsPrognostic2(double* solution,int analysis_type,int sub_analysis_type); 139 void UpdateInputsBalancedthickness( double* solution,int analysis_type,int sub_analysis_type); 140 void UpdateInputsBalancedthickness2( double* solution,int analysis_type,int sub_analysis_type); 141 void UpdateInputsBalancedvelocities( double* solution,int analysis_type,int sub_analysis_type); 142 138 143 /*}}}*/ 139 144 -
issm/trunk/src/c/objects/Sing.cpp
r3612 r3620 11 11 #include "stdio.h" 12 12 #include "./Sing.h" 13 #include "./SingVertexInput.h" 13 14 #include <string.h> 14 15 #include "../EnumDefinitions/EnumDefinitions.h" 15 16 #include "../shared/shared.h" 16 17 #include "../DataSet/DataSet.h" 18 #include "../DataSet/Inputs.h" 17 19 #include "../include/typedefs.h" 18 20 #include "../include/macros.h" … … 21 23 /*FUNCTION Sing::constructor {{{1*/ 22 24 Sing::Sing(){ 25 this->inputs=NULL; 23 26 return; 24 27 } 25 28 /*}}}*/ 26 29 /*FUNCTION Sing constructor {{{1*/ 27 Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id , ElementProperties* singproperties):30 Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id): 28 31 hnodes(sing_node_ids,1), 29 32 hmatice(&sing_matice_id,1), 30 33 hmatpar(&sing_matpar_id,1), 31 hnumpar(&sing_numpar_id,1), 32 properties(singproperties) 34 hnumpar(&sing_numpar_id,1) 33 35 { 34 36 35 37 /*all the initialization has been done by the initializer, just fill in the id: */ 36 38 this->id=sing_id; 39 this->inputs=new Inputs(); 37 40 38 41 return; … … 40 43 /*}}}*/ 41 44 /*FUNCTION Sing other constructor {{{1*/ 42 Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, ElementProperties* sing_properties):45 Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs): 43 46 hnodes(sing_hnodes), 44 47 hmatice(sing_hmatice), 45 48 hmatpar(sing_hmatpar), 46 hnumpar(sing_hnumpar), 47 properties(sing_properties) 49 hnumpar(sing_hnumpar) 48 50 { 49 51 50 52 /*all the initialization has been done by the initializer, just fill in the id: */ 51 53 this->id=sing_id; 52 54 if(sing_inputs){ 55 this->inputs=(Inputs*)sing_inputs->Copy(); 56 } 57 else{ 58 this->inputs=new Inputs(); 59 } 53 60 return; 54 61 } … … 80 87 this->hnumpar.Init(&sing_numpar_id,1); 81 88 82 / *properties: */83 sing_h=iomodel->thickness[i];84 sing_k=iomodel->drag[i]; 85 86 this->properties.Init(1,&sing_h, NULL, NULL, &sing_k, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, UNDEF,UNDEF, UNDEF,UNDEF,UNDEF);89 //intialize inputs, and add as many inputs per element as requested: 90 this->inputs=new Inputs(); 91 92 if (iomodel->thickness) this->inputs->AddInput(new SingVertexInput(ThicknessEnum,iomodel->thickness[i])); 93 if (iomodel->drag_coefficient) this->inputs->AddInput(new SingVertexInput(DragCoefficientEnum,iomodel->drag_coefficient[i])); 87 94 88 95 } … … 90 97 /*FUNCTION Sing::destructor {{{1*/ 91 98 Sing::~Sing(){ 99 delete inputs; 92 100 return; 93 101 } … … 136 144 hmatpar.DeepEcho(); 137 145 hnumpar.DeepEcho(); 138 properties.DeepEcho(); 146 printf(" inputs\n"); 147 inputs->DeepEcho(); 139 148 140 149 return; … … 160 169 hnumpar.Demarshall(&marshalled_dataset); 161 170 162 /*demarshall properties: */163 properties.Demarshall(&marshalled_dataset);171 /*demarshall inputs: */ 172 inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 164 173 165 174 /*return: */ … … 178 187 hmatpar.Echo(); 179 188 hnumpar.Echo(); 180 properties.Echo(); 181 182 return; 189 printf(" inputs\n"); 190 inputs->Echo(); 183 191 } 184 192 /*}}}*/ … … 188 196 char* marshalled_dataset=NULL; 189 197 int enum_type=0; 198 char* marshalled_inputs=NULL; 199 int marshalled_inputs_size; 190 200 191 201 /*recover marshalled_dataset: */ … … 207 217 hnumpar.Marshall(&marshalled_dataset); 208 218 209 /*Marshall properties: */ 210 properties.Marshall(&marshalled_dataset); 219 /*Marshall inputs: */ 220 marshalled_inputs_size=inputs->MarshallSize(); 221 marshalled_inputs=inputs->Marshall(); 222 memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); 223 marshalled_dataset+=marshalled_inputs_size; 224 225 xfree((void**)&marshalled_inputs); 211 226 212 227 *pmarshalled_dataset=marshalled_dataset; … … 222 237 +hmatpar.MarshallSize() 223 238 +hnumpar.MarshallSize() 224 + properties.MarshallSize()239 +inputs->MarshallSize() 225 240 +sizeof(int); //sizeof(int) for enum type 226 241 } … … 243 258 void Sing::ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type){ 244 259 245 int i; 246 const int numgrids=1; 247 int doflist[numgrids]; 248 double pressure[numgrids]; 260 int dof; 261 double pressure; 262 double thickness; 249 263 double rho_ice,g; 250 264 251 265 /*dynamic objects pointed to by hooks: */ 252 Node** nodes=NULL;253 266 Matpar* matpar=NULL; 254 267 Matice* matice=NULL; 255 Numpar* numpar=NULL;256 268 257 269 /*Get dof list on which we will plug the pressure values: */ 258 GetDofList1(&dof list[0]);270 GetDofList1(&dof); 259 271 260 272 /*recover objects from hooks: */ 261 nodes=(Node**)hnodes.deliverp();262 273 matpar=(Matpar*)hmatpar.delivers(); 263 274 matice=(Matice*)hmatice.delivers(); 264 numpar=(Numpar*)hnumpar.delivers();265 275 266 276 /*pressure is lithostatic: */ 267 277 rho_ice=matpar->GetRhoIce(); 268 278 g=matpar->GetG(); 269 pressure[0]=rho_ice*g*this->properties.h[0]; 279 inputs->GetParameterValue(&thickness,ThicknessEnum); 280 pressure=rho_ice*g*thickness; 270 281 271 282 /*plug local pressure values into global pressure vector: */ 272 VecSetValue s(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES);283 VecSetValue(p_g,dof,pressure,INSERT_VALUES); 273 284 274 285 } … … 282 293 /*}}}*/ 283 294 /*FUNCTION Sing::CostFunction {{{1*/ 284 double Sing::CostFunction( void*,int,int){295 double Sing::CostFunction( int,int){ 285 296 ISSMERROR(" not supported yet!"); 286 297 } … … 293 304 if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){ 294 305 295 CreateKMatrixDiagnosticHutter( Kgg, inputs,analysis_type,sub_analysis_type);306 CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type); 296 307 297 308 } … … 313 324 int numberofdofspernode; 314 325 315 ParameterInputs* inputs=NULL;316 317 /*recover pointers: */318 inputs=(ParameterInputs*)vinputs;319 320 326 GetDofList(&doflist[0],&numberofdofspernode); 321 327 … … 330 336 if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){ 331 337 332 CreatePVectorDiagnosticHutter( pg, inputs,analysis_type,sub_analysis_type);338 CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type); 333 339 334 340 } … … 365 371 Numpar* numpar=NULL; 366 372 367 ParameterInputs* inputs=NULL;368 369 /*recover pointers: */370 inputs=(ParameterInputs*)vinputs;371 372 373 /*recover objects from hooks: */ 373 374 node=(Node**)hnodes.deliverp(); … … 376 377 numpar=(Numpar*)hnumpar.delivers(); 377 378 378 found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,numgrids,(void**)&node[0]); 379 if(!found)ISSMERROR(" surfaceslopex missing from inputs!"); 380 found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,numgrids,(void**)&node[0]); 381 if(!found)ISSMERROR(" surfaceslopey missing from inputs!"); 379 inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum); 380 inputs->GetParameterValue(&slope[1],SurfaceSlopeyEnum); 382 381 383 382 GetDofList(&doflist[0],&numberofdofspernode); … … 391 390 n=matice->GetN(); 392 391 B=matice->GetB(); 393 thickness=this->properties.h[0];392 inputs->GetParameterValue(&thickness,ThicknessEnum); 394 393 395 394 ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0]; … … 408 407 /*}}}*/ 409 408 /*FUNCTION Sing::Du {{{1*/ 410 void Sing::Du( _p_Vec*,void*,int,int){409 void Sing::Du(Vec,int,int){ 411 410 ISSMERROR(" not supported yet!"); 412 411 } … … 504 503 /*}}}*/ 505 504 /*FUNCTION Sing::GetShelf {{{1*/ 506 intSing::GetShelf(){505 bool Sing::GetShelf(){ 507 506 ISSMERROR(" not supported yet!"); 508 507 } … … 514 513 /*}}}*/ 515 514 /*FUNCTION Sing::Gradj {{{1*/ 516 void Sing::Gradj( _p_Vec*, void*,int, int ,char*){515 void Sing::Gradj(Vec, int, int ,char*){ 517 516 ISSMERROR(" not supported yet!"); 518 517 } 519 518 /*}}}*/ 520 519 /*FUNCTION Sing::GradB {{{1*/ 521 void Sing::GradjB( _p_Vec*, void*,int,int){520 void Sing::GradjB(Vec, int,int){ 522 521 ISSMERROR(" not supported yet!"); 523 522 } 524 523 /*}}}*/ 525 524 /*FUNCTION Sing::GradjDrag {{{1*/ 526 void Sing::GradjDrag( _p_Vec*, void*,int,int){525 void Sing::GradjDrag(Vec, int,int){ 527 526 ISSMERROR(" not supported yet!"); 528 527 } … … 534 533 /*}}}*/ 535 534 /*FUNCTION Sing::Misfit {{{1*/ 536 double Sing::Misfit( void*,int,int){535 double Sing::Misfit( int,int){ 537 536 ISSMERROR(" not supported yet!"); 538 537 } … … 545 544 /*}}}*/ 546 545 /*FUNCTION Sing::SurfaceArea {{{1*/ 547 double Sing::SurfaceArea( void*,int,int){546 double Sing::SurfaceArea( int,int){ 548 547 ISSMERROR(" not supported yet!"); 549 548 } -
issm/trunk/src/c/objects/Sing.h
r3612 r3620 12 12 #include "../ModelProcessorx/IoModel.h" 13 13 #include "./Hook.h" 14 15 class Hook; 14 #include "../DataSet/Inputs.h" 16 15 17 16 class Sing: public Element{ … … 26 25 Hook hmatpar; //hook to 1 matpar 27 26 Hook hnumpar; //hook to 1 numpar 28 29 27 Inputs* inputs; 30 28 … … 39 37 /*}}}*/ 40 38 /*object management: {{{1*/ 41 void Configure( void* loads,void* nodes,void* materials,void* parameters);39 void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); 42 40 Object* copy(); 43 41 void DeepEcho(); … … 72 70 void GetBedList(double*); 73 71 void GetThicknessList(double* thickness_list); 74 void Du( _p_Vec*,void*,int,int);75 void Gradj( _p_Vec*, void*, int, int,char*);76 void GradjDrag( _p_Vec*, void*, int,int);77 void GradjB( _p_Vec*, void*,int,int);78 double Misfit( void*,int,int);79 double SurfaceArea( void*,int,int);80 double CostFunction( void*,int,int);72 void Du(Vec,int,int); 73 void Gradj(Vec, int, int,char*); 74 void GradjDrag(Vec , int,int); 75 void GradjB(Vec, int,int); 76 double Misfit(int,int); 77 double SurfaceArea(int,int); 78 double CostFunction(int,int); 81 79 double MassFlux(double* segment,double* ug); 82 80 /*}}}*/ -
issm/trunk/src/c/objects/Tria.cpp
r3612 r3620 30 30 } 31 31 /*}}}*/ 32 /*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id , int numpar_id){{{1*/33 Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id , int tria_numpar_id):32 /*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id){{{1*/ 33 Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id): 34 34 hnodes(tria_node_ids,3), 35 35 hmatice(&tria_matice_id,1), 36 hmatpar(&tria_matpar_id,1), 37 hnumpar(&tria_numpar_id,1) 36 hmatpar(&tria_matpar_id,1) 38 37 { 39 38 40 39 /*all the initialization has been done by the initializer, just fill in the id: */ 41 40 this->id=tria_id; 41 this->parameters=NULL; 42 42 this->inputs=new Inputs(); 43 43 44 return; 45 } 46 /*}}}*/ 47 /*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, Inputs* tria_inputs) {{{1*/ 48 Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs): 44 } 45 /*}}}*/ 46 /*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* tria_inputs) {{{1*/ 47 Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* tria_parameters, Inputs* tria_inputs): 49 48 hnodes(tria_hnodes), 50 49 hmatice(tria_hmatice), 51 hmatpar(tria_hmatpar), 52 hnumpar(tria_hnumpar) 50 hmatpar(tria_hmatpar) 53 51 { 54 52 … … 61 59 this->inputs=new Inputs(); 62 60 } 63 64 return;61 /*point parameters: */ 62 this->parameters=tria_parameters; 65 63 } 66 64 /*}}}*/ … … 72 70 int tria_matice_id; 73 71 int tria_matpar_id; 74 int tria_numpar_id;75 72 double nodeinputs[3]; 76 73 … … 94 91 tria_matice_id=index+1; //refers to the corresponding ice material object 95 92 tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object 96 tria_numpar_id=1; //refers to numerical parameters object97 93 98 94 this->hnodes.Init(tria_node_ids,3); 99 95 this->hmatice.Init(&tria_matice_id,1); 100 96 this->hmatpar.Init(&tria_matpar_id,1); 101 this->hnumpar.Init(&tria_numpar_id,1);102 97 103 98 //intialize inputs, and add as many inputs per element as requested: … … 143 138 if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); 144 139 140 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. 141 this->parameters=NULL; 142 145 143 146 144 } … … 149 147 Tria::~Tria(){ 150 148 delete inputs; 149 this->parameters=NULL; 151 150 return; 152 151 } … … 155 154 /*Object management: */ 156 155 /*FUNCTION Tria::Configure {{{1*/ 157 void Tria::Configure( void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){156 void Tria::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ 158 157 159 158 int i; 160 161 DataSet* loadsin=NULL;162 DataSet* nodesin=NULL;163 DataSet* materialsin=NULL;164 DataSet* parametersin=NULL;165 166 /*Recover pointers :*/167 loadsin=(DataSet*)ploadsin;168 nodesin=(DataSet*)pnodesin;169 materialsin=(DataSet*)pmaterialsin;170 parametersin=(DataSet*)pparametersin;171 159 172 160 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective … … 175 163 hmatice.configure(materialsin); 176 164 hmatpar.configure(materialsin); 177 hnumpar.configure(parametersin); 165 166 /*point parameters to real dataset: */ 167 this->parameters=parametersin; 178 168 179 169 } … … 182 172 Object* Tria::copy() { 183 173 184 return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar, &this->hnumpar,this->inputs);174 return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs); 185 175 186 176 } … … 205 195 hmatice.Demarshall(&marshalled_dataset); 206 196 hmatpar.Demarshall(&marshalled_dataset); 207 hnumpar.Demarshall(&marshalled_dataset);208 197 209 198 /*demarshall inputs: */ 210 199 inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 200 201 /*parameters: may not exist even yet, so let Configure handle it: */ 202 this->parameters=NULL; 211 203 212 204 /*return: */ … … 224 216 hmatice.DeepEcho(); 225 217 hmatpar.DeepEcho(); 226 hnumpar.DeepEcho(); 218 printf(" parameters\n"); 219 parameters->DeepEcho(); 227 220 printf(" inputs\n"); 228 221 inputs->DeepEcho(); 229 222 230 223 return; 231 224 } … … 240 233 hmatice.Echo(); 241 234 hmatpar.Echo(); 242 hnumpar.Echo(); 235 printf(" parameters\n"); 236 parameters->Echo(); 243 237 printf(" inputs\n"); 244 238 inputs->Echo(); 245 246 return;247 239 } 248 240 /*}}}*/ … … 271 263 hmatice.Marshall(&marshalled_dataset); 272 264 hmatpar.Marshall(&marshalled_dataset); 273 hnumpar.Marshall(&marshalled_dataset);274 265 275 266 /*Marshall inputs: */ … … 279 270 marshalled_dataset+=marshalled_inputs_size; 280 271 272 /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/ 273 274 xfree((void**)&marshalled_inputs); 275 281 276 *pmarshalled_dataset=marshalled_dataset; 282 277 return; … … 290 285 +hmatice.MarshallSize() 291 286 +hmatpar.MarshallSize() 292 +hnumpar.MarshallSize()293 287 +inputs->MarshallSize() 294 288 +sizeof(int); //sizeof(int) for enum type … … 312 306 Matpar* matpar=NULL; 313 307 Matice* matice=NULL; 314 Numpar* numpar=NULL;315 308 316 309 /*recover objects from hooks: */ … … 318 311 matpar=(Matpar*)hmatpar.delivers(); 319 312 matice=(Matice*)hmatice.delivers(); 320 numpar=(Numpar*)hnumpar.delivers();321 313 322 314 /*Update internal data if inputs holds new values: */ … … 468 460 Matpar* matpar=NULL; 469 461 Matice* matice=NULL; 470 Numpar* numpar=NULL;471 462 472 463 /*recover objects from hooks: */ … … 474 465 matpar=(Matpar*)hmatpar.delivers(); 475 466 matice=(Matice*)hmatice.delivers(); 476 numpar=(Numpar*)hnumpar.delivers();477 467 478 468 /*plug local pressure values into global pressure vector: */ … … 497 487 Matpar* matpar=NULL; 498 488 Matice* matice=NULL; 499 Numpar* numpar=NULL;500 489 501 490 /*recover objects from hooks: */ … … 503 492 matpar=(Matpar*)hmatpar.delivers(); 504 493 matice=(Matice*)hmatice.delivers(); 505 numpar=(Numpar*)hnumpar.delivers();506 494 507 495 /*Get dof list on which we will plug the pressure values: */ … … 536 524 Matpar* matpar=NULL; 537 525 Matice* matice=NULL; 538 Numpar* numpar=NULL;539 526 540 527 /*recover objects from hooks: */ … … 542 529 matpar=(Matpar*)hmatpar.delivers(); 543 530 matice=(Matice*)hmatice.delivers(); 544 numpar=(Numpar*)hnumpar.delivers();545 531 546 532 /*plug local pressure values into global pressure vector: */ … … 583 569 double dk[NDOF2]; 584 570 double dB[NDOF2]; 571 char* control_type=NULL; 572 double cm_noisedmp; 585 573 586 574 /* Jacobian: */ … … 591 579 Matpar* matpar=NULL; 592 580 Matice* matice=NULL; 593 Numpar* numpar=NULL;594 581 595 582 /*inputs: */ … … 608 595 matpar=(Matpar*)hmatpar.delivers(); 609 596 matice=(Matice*)hmatice.delivers(); 610 numpar=(Numpar*)hnumpar.delivers(); 597 598 /*retrieve some parameters: */ 599 this->parameters->FindParam(&control_type,"control_type"); 600 this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); 611 601 612 602 /* Get node coordinates and dof list: */ … … 631 621 632 622 /*Add Tikhonov regularization term to misfit*/ 633 if (strcmp( numpar->control_type,"drag")==0){623 if (strcmp(control_type,"drag")==0){ 634 624 if (shelf){ 635 625 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); 636 Jelem+= numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight;626 Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; 637 627 638 628 } 639 629 } 640 else if (strcmp( numpar->control_type,"B")==0){630 else if (strcmp(control_type,"B")==0){ 641 631 inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBEnum); 642 Jelem+= numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;632 Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; 643 633 } 644 634 else{ 645 ISSMERROR("%s%s","unsupported control type: ", numpar->control_type);635 ISSMERROR("%s%s","unsupported control type: ",control_type); 646 636 } 647 637 … … 652 642 xfree((void**)&third_gauss_area_coord); 653 643 xfree((void**)&gauss_weights); 644 xfree((void**)&control_type); 654 645 655 646 /*Return: */ … … 763 754 Matpar* matpar=NULL; 764 755 Matice* matice=NULL; 765 Numpar* numpar=NULL; 756 757 /*parameters: */ 758 double artdiff; 766 759 767 760 … … 770 763 matpar=(Matpar*)hmatpar.delivers(); 771 764 matice=(Matice*)hmatice.delivers(); 772 numpar=(Numpar*)hnumpar.delivers();773 765 774 766 /* Get node coordinates and dof list: */ … … 780 772 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); 781 773 774 /*retrieve some parameters: */ 775 this->parameters->FindParam(&artdiff,"artdiff"); 776 782 777 //Create Artificial diffusivity once for all if requested 783 if( numpar->artdiff){778 if(artdiff){ 784 779 //Get the Jacobian determinant 785 780 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; … … 848 843 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j]; 849 844 850 if( numpar->artdiff){845 if(artdiff){ 851 846 852 847 /* Compute artificial diffusivity */ … … 918 913 Matpar* matpar=NULL; 919 914 Matice* matice=NULL; 920 Numpar* numpar=NULL;921 915 922 916 … … 925 919 matpar=(Matpar*)hmatpar.delivers(); 926 920 matice=(Matice*)hmatice.delivers(); 927 numpar=(Numpar*)hnumpar.delivers();928 921 929 922 /* Get node coordinates and dof list: */ … … 1034 1027 int found=0; 1035 1028 1029 /*parameters: */ 1030 double artdiff; 1031 1036 1032 /*dynamic objects pointed to by hooks: */ 1037 1033 Node** nodes=NULL; 1038 1034 Matpar* matpar=NULL; 1039 1035 Matice* matice=NULL; 1040 Numpar* numpar=NULL;1041 1036 1042 1037 /*recover objects from hooks: */ … … 1044 1039 matpar=(Matpar*)hmatpar.delivers(); 1045 1040 matice=(Matice*)hmatice.delivers(); 1046 numpar=(Numpar*)hnumpar.delivers();1047 1041 1048 1042 /*Recover velocity: */ 1049 1043 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); 1050 1044 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); 1045 1046 /*retrieve some parameters: */ 1047 this->parameters->FindParam(&artdiff,"artdiff"); 1051 1048 1052 1049 /* Get node coordinates and dof list: */ … … 1075 1072 1076 1073 //Create Artificial diffusivity once for all if requested 1077 if( numpar->artdiff){1074 if(artdiff){ 1078 1075 //Get the Jacobian determinant 1079 1076 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; … … 1131 1128 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j]; 1132 1129 1133 if( numpar->artdiff){1130 if(artdiff){ 1134 1131 1135 1132 /* Compute artificial diffusivity */ … … 1195 1192 double B[3][numdof]; 1196 1193 double Bprime[3][numdof]; 1197 double D[3][3]={ { 0,0,0 },{0,0,0},{0,0,0}};// material matrix, simple scalar matrix.1194 double D[3][3]={0.0}; // material matrix, simple scalar matrix. 1198 1195 double D_scalar; 1199 1196 1197 /*parameters: */ 1198 double viscosity_overshoot; 1199 1200 1200 /* local element matrices: */ 1201 double Ke_gg[numdof][numdof] ; //local element stiffness matrix1201 double Ke_gg[numdof][numdof]={0.0}; 1202 1202 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. 1203 1203 … … 1212 1212 Matpar* matpar=NULL; 1213 1213 Matice* matice=NULL; 1214 Numpar* numpar=NULL;1215 1214 1216 1215 /*inputs: */ … … 1220 1219 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1221 1220 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 1221 1222 /*retrieve some parameters: */ 1223 this->parameters->FindParam(&viscosity_overshoot,"viscosity_overshoot"); 1222 1224 1223 1225 /*First, if we are on water, return empty matrix: */ … … 1228 1230 matpar=(Matpar*)hmatpar.delivers(); 1229 1231 matice=(Matice*)hmatice.delivers(); 1230 numpar=(Numpar*)hnumpar.delivers();1231 1232 1232 1233 /* Get node coordinates and dof list: */ 1233 1234 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 1234 1235 GetDofList(&doflist[0],&numberofdofspernode); 1235 1236 /* Set Ke_gg to 0: */1237 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;1238 1236 1239 1237 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 1265 1263 /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 1266 1264 onto this scalar matrix, so that we win some computational time: */ 1267 newviscosity=viscosity+ numpar->viscosity_overshoot*(viscosity-oldviscosity);1265 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 1268 1266 D_scalar=newviscosity*thickness*gauss_weight*Jdet; 1269 1267 … … 1333 1331 1334 1332 /* local element matrices: */ 1335 double Ke_gg[numdof][numdof] ; //local element stiffness matrix1333 double Ke_gg[numdof][numdof]={0.0}; 1336 1334 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag 1337 1335 … … 1365 1363 Matpar* matpar=NULL; 1366 1364 Matice* matice=NULL; 1367 Numpar* numpar=NULL;1368 1365 1369 1366 … … 1372 1369 matpar=(Matpar*)hmatpar.delivers(); 1373 1370 matice=(Matice*)hmatice.delivers(); 1374 numpar=(Numpar*)hnumpar.delivers();1375 1371 1376 1372 /*retrieve inputs :*/ … … 1381 1377 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 1382 1378 GetDofList(&doflist[0],&numberofdofspernode); 1383 1384 /* Set Ke_gg to 0: */1385 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;1386 1379 1387 1380 if (shelf){ … … 1522 1515 1523 1516 /* local element matrices: */ 1524 double Ke_gg[numdof][numdof] ; //local element stiffness matrix1517 double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 1525 1518 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. 1526 1519 … … 1529 1522 Matpar* matpar=NULL; 1530 1523 Matice* matice=NULL; 1531 Numpar* numpar=NULL;1532 1524 1533 1525 /*recover objects from hooks: */ … … 1535 1527 matpar=(Matpar*)hmatpar.delivers(); 1536 1528 matice=(Matice*)hmatice.delivers(); 1537 numpar=(Numpar*)hnumpar.delivers();1538 1529 1539 1530 /* Get node coordinates and dof list: */ 1540 1531 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 1541 1532 GetDofList(&doflist[0],&numberofdofspernode); 1542 1543 /* Set Ke_gg to 0: */1544 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;1545 1533 1546 1534 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 1655 1643 Matpar* matpar=NULL; 1656 1644 Matice* matice=NULL; 1657 Numpar* numpar=NULL;1658 1645 1659 1646 /*recover objects from hooks: */ … … 1661 1648 matpar=(Matpar*)hmatpar.delivers(); 1662 1649 matice=(Matice*)hmatice.delivers(); 1663 numpar=(Numpar*)hnumpar.delivers();1664 1650 1665 1651 /*Recover constants of ice */ … … 1759 1745 double K[2][2]={0.0}; 1760 1746 double KDL[2][2]={0.0}; 1761 double dt;1762 1747 int dofs[2]={0,1}; 1763 1748 int found; 1749 1750 /*parameters: */ 1751 double dt; 1752 double artdiff; 1764 1753 1765 1754 /*dynamic objects pointed to by hooks: */ … … 1767 1756 Matpar* matpar=NULL; 1768 1757 Matice* matice=NULL; 1769 Numpar* numpar=NULL;1770 1758 1771 1759 /*recover objects from hooks: */ … … 1773 1761 matpar=(Matpar*)hmatpar.delivers(); 1774 1762 matice=(Matice*)hmatice.delivers(); 1775 numpar=(Numpar*)hnumpar.delivers(); 1776 1777 /*recover dt: */1778 dt=numpar->dt;1763 1764 /*retrieve some parameters: */ 1765 this->parameters->FindParam(&dt,"dt"); 1766 this->parameters->FindParam(&artdiff,"artdiff"); 1779 1767 1780 1768 /* Get node coordinates and dof list: */ … … 1787 1775 1788 1776 //Create Artificial diffusivity once for all if requested 1789 if( numpar->artdiff==1){1777 if(artdiff==1){ 1790 1778 //Get the Jacobian determinant 1791 1779 gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; … … 1867 1855 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j]; 1868 1856 1869 if( numpar->artdiff==1){1857 if(artdiff==1){ 1870 1858 1871 1859 /* Compute artificial diffusivity */ … … 1934 1922 /*input parameters for structural analysis (diagnostic): */ 1935 1923 double vx,vy; 1936 double dt;1937 1924 int dofs[1]={0}; 1938 1925 int found; 1926 1927 /*parameters: */ 1928 double dt; 1939 1929 1940 1930 /*dynamic objects pointed to by hooks: */ … … 1942 1932 Matpar* matpar=NULL; 1943 1933 Matice* matice=NULL; 1944 Numpar* numpar=NULL;1945 1934 1946 1935 /*recover objects from hooks: */ … … 1948 1937 matpar=(Matpar*)hmatpar.delivers(); 1949 1938 matice=(Matice*)hmatice.delivers(); 1950 numpar=(Numpar*)hnumpar.delivers(); 1951 1952 /*recover dt: */ 1953 dt=numpar->dt; 1939 1940 /*retrieve some parameters: */ 1941 this->parameters->FindParam(&dt,"dt"); 1954 1942 1955 1943 /* Get node coordinates and dof list: */ … … 2050 2038 2051 2039 /* local element matrices: */ 2052 double Ke_gg[numdof][numdof] ; //local element stiffness matrix2040 double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 2053 2041 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. 2054 2042 … … 2059 2047 Matpar* matpar=NULL; 2060 2048 Matice* matice=NULL; 2061 Numpar* numpar=NULL;2062 2049 2063 2050 /*recover objects from hooks: */ … … 2065 2052 matpar=(Matpar*)hmatpar.delivers(); 2066 2053 matice=(Matice*)hmatice.delivers(); 2067 numpar=(Numpar*)hnumpar.delivers();2068 2054 2069 2055 /* Get node coordinates and dof list: */ 2070 2056 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2071 2057 GetDofList(&doflist[0],&numberofdofspernode); 2072 2073 /* Set Ke_gg to 0: */2074 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;2075 2058 2076 2059 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 2133 2116 double rho_ice; 2134 2117 double heatcapacity; 2135 double dt;2136 2118 2137 2119 int num_gauss,ig; … … 2151 2133 double D_scalar; 2152 2134 2135 /*parameters: */ 2136 double dt; 2137 2153 2138 /*dynamic objects pointed to by hooks: */ 2154 2139 Node** nodes=NULL; 2155 2140 Matpar* matpar=NULL; 2156 2141 Matice* matice=NULL; 2157 Numpar* numpar=NULL;2158 2142 2159 2143 /*recover objects from hooks: */ … … 2161 2145 matpar=(Matpar*)hmatpar.delivers(); 2162 2146 matice=(Matice*)hmatice.delivers(); 2163 numpar=(Numpar*)hnumpar.delivers(); 2164 2165 /*recover extra inputs from users, dt: */ 2166 dt=numpar->dt; 2147 2148 /*retrieve some parameters: */ 2149 this->parameters->FindParam(&dt,"dt"); 2167 2150 2168 2151 /* Get node coordinates and dof list: */ … … 2304 2287 Matpar* matpar=NULL; 2305 2288 Matice* matice=NULL; 2306 Numpar* numpar=NULL;2307 2289 2308 2290 /*recover objects from hooks: */ … … 2310 2292 matpar=(Matpar*)hmatpar.delivers(); 2311 2293 matice=(Matice*)hmatice.delivers(); 2312 numpar=(Numpar*)hnumpar.delivers();2313 2294 2314 2295 /* Get node coordinates and dof list: */ … … 2391 2372 Matpar* matpar=NULL; 2392 2373 Matice* matice=NULL; 2393 Numpar* numpar=NULL;2394 2374 2395 2375 /*recover objects from hooks: */ … … 2397 2377 matpar=(Matpar*)hmatpar.delivers(); 2398 2378 matice=(Matice*)hmatice.delivers(); 2399 numpar=(Numpar*)hnumpar.delivers();2400 2379 2401 2380 /* Get node coordinates and dof list: */ … … 2478 2457 Matpar* matpar=NULL; 2479 2458 Matice* matice=NULL; 2480 Numpar* numpar=NULL;2481 2459 2482 2460 /*recover objects from hooks: */ … … 2484 2462 matpar=(Matpar*)hmatpar.delivers(); 2485 2463 matice=(Matice*)hmatice.delivers(); 2486 numpar=(Numpar*)hnumpar.delivers();2487 2488 /*recover objects from hooks: */2489 nodes=(Node**)hnodes.deliverp();2490 matpar=(Matpar*)hmatpar.delivers();2491 matice=(Matice*)hmatice.delivers();2492 numpar=(Numpar*)hnumpar.delivers();2493 2464 2494 2465 /* Get node coordinates and dof list: */ … … 2562 2533 2563 2534 /*element vector at the gaussian points: */ 2564 double pe_g[numdof] ;2535 double pe_g[numdof]={0.0}; 2565 2536 double pe_g_gaussian[numdof]; 2566 2537 … … 2578 2549 Matpar* matpar=NULL; 2579 2550 Matice* matice=NULL; 2580 Numpar* numpar=NULL;2581 2551 2582 2552 /*recover objects from hooks: */ … … 2584 2554 matpar=(Matpar*)hmatpar.delivers(); 2585 2555 matice=(Matice*)hmatice.delivers(); 2586 numpar=(Numpar*)hnumpar.delivers();2587 2556 2588 2557 /* Get node coordinates and dof list: */ 2589 2558 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2590 2559 GetDofList(&doflist[0],&numberofdofspernode); 2591 2592 /* Set pe_g to 0: */2593 for(i=0;i<numdof;i++) pe_g[i]=0.0;2594 2560 2595 2561 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 2680 2646 2681 2647 /*element vector at the gaussian points: */ 2682 double pe_g[numdof] ;2648 double pe_g[numdof]={0.0}; 2683 2649 double pe_g_gaussian[numdof]; 2684 2650 … … 2690 2656 Matpar* matpar=NULL; 2691 2657 Matice* matice=NULL; 2692 Numpar* numpar=NULL;2693 2658 2694 2659 /*inputs: */ … … 2707 2672 matpar=(Matpar*)hmatpar.delivers(); 2708 2673 matice=(Matice*)hmatice.delivers(); 2709 numpar=(Numpar*)hnumpar.delivers();2710 2674 2711 2675 /* Get node coordinates and dof list: */ 2712 2676 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2713 2677 GetDofList(&doflist[0],&numberofdofspernode); 2714 2715 /* Set pe_g to 0: */2716 for(i=0;i<numdof;i++) pe_g[i]=0.0;2717 2678 2718 2679 … … 2812 2773 double melting_g; 2813 2774 double thickness_g; 2775 2776 /*parameters: */ 2814 2777 double dt; 2815 2778 … … 2818 2781 Matpar* matpar=NULL; 2819 2782 Matice* matice=NULL; 2820 Numpar* numpar=NULL;2821 2783 2822 2784 /*recover objects from hooks: */ … … 2824 2786 matpar=(Matpar*)hmatpar.delivers(); 2825 2787 matice=(Matice*)hmatice.delivers(); 2826 numpar=(Numpar*)hnumpar.delivers(); 2827 2828 /*recover dt: */ 2829 dt=numpar->dt; 2788 2789 /*retrieve some parameters: */ 2790 this->parameters->FindParam(&dt,"dt"); 2830 2791 2831 2792 /* Get node coordinates and dof list: */ … … 2909 2870 Matpar* matpar=NULL; 2910 2871 Matice* matice=NULL; 2911 Numpar* numpar=NULL;2912 2872 2913 2873 /*recover objects from hooks: */ … … 2915 2875 matpar=(Matpar*)hmatpar.delivers(); 2916 2876 matice=(Matice*)hmatice.delivers(); 2917 numpar=(Numpar*)hnumpar.delivers(); 2918 2919 /*recover dt: */ 2920 dt=numpar->dt; 2921 2877 2878 /*retrieve some parameters: */ 2879 this->parameters->FindParam(&dt,"dt"); 2922 2880 2923 2881 /* Get node coordinates and dof list: */ … … 2993 2951 2994 2952 /*element vector at the gaussian points: */ 2995 double pe_g[numdof] ;2953 double pe_g[numdof]={0.0}; 2996 2954 double pe_g_gaussian[numdof]; 2997 2955 double slope[2]; … … 3001 2959 Matpar* matpar=NULL; 3002 2960 Matice* matice=NULL; 3003 Numpar* numpar=NULL;3004 2961 3005 2962 /*recover objects from hooks: */ … … 3007 2964 matpar=(Matpar*)hmatpar.delivers(); 3008 2965 matice=(Matice*)hmatice.delivers(); 3009 numpar=(Numpar*)hnumpar.delivers();3010 2966 3011 2967 /* Get node coordinates and dof list: */ 3012 2968 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 3013 2969 GetDofList(&doflist[0],&numberofdofspernode); 3014 3015 /* Set pe_g to 0: */3016 for(i=0;i<numdof;i++) pe_g[i]=0.0;3017 2970 3018 2971 … … 3112 3065 Matpar* matpar=NULL; 3113 3066 Matice* matice=NULL; 3114 Numpar* numpar=NULL;3115 3067 3116 3068 /*recover objects from hooks: */ … … 3118 3070 matpar=(Matpar*)hmatpar.delivers(); 3119 3071 matice=(Matice*)hmatice.delivers(); 3120 numpar=(Numpar*)hnumpar.delivers();3121 3072 3122 3073 /* Get node coordinates and dof list: */ … … 3132 3083 beta=matpar->GetBeta(); 3133 3084 meltingpoint=matpar->GetMeltingPoint(); 3134 dt=numpar->dt; 3135 3085 3086 /*retrieve some solution parameters: */ 3087 this->parameters->FindParam(&dt,"dt"); 3136 3088 3137 3089 /* Ice/ocean heat exchange flux on ice shelf base */ … … 3231 3183 Matpar* matpar=NULL; 3232 3184 Matice* matice=NULL; 3233 Numpar* numpar=NULL;3234 3185 3235 3186 /*recover objects from hooks: */ … … 3237 3188 matpar=(Matpar*)hmatpar.delivers(); 3238 3189 matice=(Matice*)hmatice.delivers(); 3239 numpar=(Numpar*)hnumpar.delivers();3240 3190 3241 3191 /*retrieve inputs :*/ … … 3249 3199 rho_ice=matpar->GetRhoIce(); 3250 3200 heatcapacity=matpar->GetHeatCapacity(); 3251 dt=numpar->dt; 3201 3202 /*retrieve some parameters: */ 3203 this->parameters->FindParam(&dt,"dt"); 3252 3204 3253 3205 … … 3369 3321 double obs_velocity_mag,velocity_mag; 3370 3322 double dux,duy; 3323 double meanvel, epsvel; 3371 3324 3372 3325 /*element vector : */ 3373 double due_g[numdof]={0 ,0,0,0,0,0};3326 double due_g[numdof]={0.0}; 3374 3327 double due_g_gaussian[numdof]; 3375 3328 … … 3391 3344 Matpar* matpar=NULL; 3392 3345 Matice* matice=NULL; 3393 Numpar* numpar=NULL;3394 3346 3395 3347 /*recover objects from hooks: */ … … 3397 3349 matpar=(Matpar*)hmatpar.delivers(); 3398 3350 matice=(Matice*)hmatice.delivers(); 3399 numpar=(Numpar*)hnumpar.delivers(); 3351 3352 /*retrieve some parameters: */ 3353 this->parameters->FindParam(&meanvel,"meanvel"); 3354 this->parameters->FindParam(&epsvel,"epsvel"); 3400 3355 3401 3356 /* Get node coordinates and dof list: */ … … 3460 3415 */ 3461 3416 for (i=0;i<numgrids;i++){ 3462 scalex=pow( numpar->meanvel/(obs_vx_list[i]+numpar->epsvel),2);3463 scaley=pow( numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),2);3417 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2); 3418 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2); 3464 3419 if(obs_vx_list[i]==0)scalex=0; 3465 3420 if(obs_vy_list[i]==0)scaley=0; … … 3482 3437 */ 3483 3438 for (i=0;i<numgrids;i++){ 3484 velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+ numpar->epsvel; //epsvel to avoid velocity being nil.3485 obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+ numpar->epsvel; //epsvel to avoid observed velocity being nil.3486 scale=-8*pow( numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);3439 velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil. 3440 obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil. 3441 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 3487 3442 dux_list[i]=scale*vx_list[i]; 3488 3443 duy_list[i]=scale*vy_list[i]; … … 3501 3456 */ 3502 3457 for (i=0;i<numgrids;i++){ 3503 scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+ numpar->epsvel);3458 scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel); 3504 3459 dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]); 3505 3460 duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]); … … 3518 3473 */ 3519 3474 for (i=0;i<numgrids;i++){ 3520 dux_list[i] = - pow( numpar->meanvel,(double)2)*(3521 log((fabs(vx_list[i])+ numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)) * 1/(vx_list[i]+numpar->epsvel));3522 duy_list[i] = - pow( numpar->meanvel,(double)2)*(3523 log((fabs(vy_list[i])+ numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)) * 1/(vy_list[i]+numpar->epsvel));3475 dux_list[i] = - pow(meanvel,(double)2)*( 3476 log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel)); 3477 duy_list[i] = - pow(meanvel,(double)2)*( 3478 log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel)); 3524 3479 } 3525 3480 } … … 4228 4183 4229 4184 /*element vector at the gaussian points: */ 4230 double grade_g[numgrids] ;4185 double grade_g[numgrids]={0.0}; 4231 4186 double grade_g_gaussian[numgrids]; 4232 4187 … … 4252 4207 double dB[NDOF2]; 4253 4208 double B_gauss; 4209 4210 /*parameters: */ 4211 double cm_noisedmp; 4212 double cm_mindmp_slope; 4213 double cm_mindmp_value; 4214 double cm_maxdmp_value; 4215 double cm_maxdmp_slope; 4254 4216 4255 4217 /*dynamic objects pointed to by hooks: */ … … 4257 4219 Matpar* matpar=NULL; 4258 4220 Matice* matice=NULL; 4259 Numpar* numpar=NULL;4260 4221 4261 4222 /*recover objects from hooks: */ … … 4263 4224 matpar=(Matpar*)hmatpar.delivers(); 4264 4225 matice=(Matice*)hmatice.delivers(); 4265 numpar=(Numpar*)hnumpar.delivers(); 4226 4227 /*retrieve some parameters: */ 4228 this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); 4229 this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value"); 4230 this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope"); 4231 this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); 4232 this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); 4266 4233 4267 4234 /* Get node coordinates and dof list: */ 4268 4235 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 4269 4236 GetDofList1(&doflist1[0]); 4270 4271 /* Set grade_g to 0: */4272 for(i=0;i<numgrids;i++) grade_g[i]=0.0;4273 4237 4274 4238 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 4317 4281 4318 4282 //Add regularization term 4319 grade_g_gaussian[i]-= numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]);4283 grade_g_gaussian[i]-=cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]); 4320 4284 4321 4285 //min dampening 4322 if(B_gauss< numpar->cm_mindmp_value){4323 grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];4286 if(B_gauss<cm_mindmp_value){ 4287 grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; 4324 4288 } 4325 4289 4326 4290 //max dampening 4327 if(B_gauss> numpar->cm_maxdmp_value){4328 grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];4291 if(B_gauss>cm_maxdmp_value){ 4292 grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; 4329 4293 } 4330 4294 … … 4394 4358 4395 4359 /*element vector at the gaussian points: */ 4396 double grade_g[numgrids]={0 ,0,0};4360 double grade_g[numgrids]={0.0}; 4397 4361 double grade_g_gaussian[numgrids]; 4398 4362 … … 4408 4372 /*inputs: */ 4409 4373 bool shelf; 4374 4375 /*parameters: */ 4376 double cm_noisedmp; 4377 double cm_mindmp_slope; 4378 double cm_mindmp_value; 4379 double cm_maxdmp_value; 4380 double cm_maxdmp_slope; 4410 4381 4411 4382 /*dynamic objects pointed to by hooks: */ … … 4413 4384 Matpar* matpar=NULL; 4414 4385 Matice* matice=NULL; 4415 Numpar* numpar=NULL;4416 4386 4417 4387 /*recover objects from hooks: */ … … 4419 4389 matpar=(Matpar*)hmatpar.delivers(); 4420 4390 matice=(Matice*)hmatice.delivers(); 4421 numpar=(Numpar*)hnumpar.delivers();4422 4391 4423 4392 /*retrieve inputs :*/ 4424 4393 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 4394 4395 /*retrieve some parameters: */ 4396 this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); 4397 this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value"); 4398 this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope"); 4399 this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); 4400 this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); 4425 4401 4426 4402 /*Get out if shelf*/ … … 4520 4496 4521 4497 //noise dampening d/dki(1/2*(dk/dx)^2) 4522 grade_g_gaussian[i]+=- numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);4498 grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); 4523 4499 4524 4500 //min dampening 4525 if(drag< numpar->cm_mindmp_value){4526 grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];4501 if(drag<cm_mindmp_value){ 4502 grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; 4527 4503 } 4528 4504 4529 4505 //max dampening 4530 if(drag> numpar->cm_maxdmp_value){4531 grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];4506 if(drag>cm_maxdmp_value){ 4507 grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; 4532 4508 } 4533 4509 } … … 4564 4540 double vy_list[numgrids]; 4565 4541 double vz_list[numgrids]; 4566 double vxvyvz_list[numgrids][3];4567 double adjx_list[numgrids];4568 double adjy_list[numgrids];4569 double adjz_list[numgrids];4570 double adjxyz_list[numgrids][3];4571 4572 4542 double drag; 4573 int dofs1[1]={0}; 4574 int dofs3[3]={0,1,2}; 4543 double thickness_list[numgrids]; 4544 double bed_list[numgrids]; 4545 double dragcoefficient_list[numgrids]; 4546 double drag_p,drag_q; 4547 double alpha_complement_list[numgrids]; 4548 double alpha_complement; 4575 4549 4576 4550 /* gaussian points: */ … … 4582 4556 double gauss_weight; 4583 4557 double gauss_l1l2l3[3]; 4558 double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; 4584 4559 4585 4560 /* parameters: */ … … 4591 4566 double dk[NDOF2]; 4592 4567 4593 /*drag: */4594 double pcoeff,qcoeff;4595 double alpha_complement_list[numgrids];4596 double alpha_complement;4597 4598 4568 /*element vector at the gaussian points: */ 4599 double grade_g[numgrids] ;4569 double grade_g[numgrids]={0.0}; 4600 4570 double grade_g_gaussian[numgrids]; 4601 4571 … … 4612 4582 bool shelf; 4613 4583 int drag_type; 4584 4585 /*parameters: */ 4586 double cm_noisedmp; 4587 double cm_mindmp_slope; 4588 double cm_mindmp_value; 4589 double cm_maxdmp_value; 4590 double cm_maxdmp_slope; 4614 4591 4615 4592 /*dynamic objects pointed to by hooks: */ … … 4617 4594 Matpar* matpar=NULL; 4618 4595 Matice* matice=NULL; 4619 Numpar* numpar=NULL;4620 4596 4621 4597 /*retrieve inputs :*/ 4622 4598 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 4623 4599 inputs->GetParameterValue(&drag_type,DragTypeEnum); 4600 4601 /*retrieve some parameters: */ 4602 this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); 4603 this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value"); 4604 this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope"); 4605 this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); 4606 this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); 4624 4607 4625 4608 /*Get out if shelf*/ … … 4630 4613 matpar=(Matpar*)hmatpar.delivers(); 4631 4614 matice=(Matice*)hmatice.delivers(); 4632 numpar=(Numpar*)hnumpar.delivers();4633 4615 4634 4616 /* Get node coordinates and dof list: */ … … 4636 4618 GetDofList1(&doflist1[0]); 4637 4619 4638 /* Set grade_g to 0: */ 4639 for(i=0;i<numgrids;i++) grade_g[i]=0.0; 4640 4641 /* recover input parameters: */ 4642 inputs->Recover("drag",&this->properties.k[0],1,dofs1,numgrids,(void**)nodes); 4643 inputs->Recover("bed",&this->properties.b[0],1,dofs1,numgrids,(void**)nodes); 4644 inputs->Recover("thickness",&this->properties.h[0],1,dofs1,numgrids,(void**)nodes); 4645 if(!inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs3,numgrids,(void**)nodes)){ 4646 ISSMERROR("missing velocity input parameter"); 4647 } 4648 if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,numgrids,(void**)nodes)){ 4649 ISSMERROR("missing adjoint input parameter"); 4650 } 4651 4652 /*Initialize parameter lists: */ 4653 for(i=0;i<numgrids;i++){ 4654 vx_list[i]=vxvyvz_list[i][0]; 4655 vy_list[i]=vxvyvz_list[i][1]; 4656 vz_list[i]=vxvyvz_list[i][2]; 4657 adjx_list[i]=adjxyz_list[i][0]; 4658 adjy_list[i]=adjxyz_list[i][1]; 4659 adjz_list[i]=adjxyz_list[i][2]; 4620 /*Build alpha_complement_list: */ 4621 if (drag_type==2){ 4622 4623 /*Allocate friction object: */ 4624 Friction* friction=NewFriction(); 4625 4626 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); 4627 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); 4628 inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],3,VzAverageEnum); 4629 inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum); 4630 inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum); 4631 inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum); 4632 inputs->GetParameterValue(&drag_p,DragPEnum); 4633 inputs->GetParameterValue(&drag_q,DragQEnum); 4634 4635 4636 /*Initialize all fields: */ 4637 friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); 4638 strcpy(friction->element_type,"2d"); 4639 4640 friction->gravity=matpar->GetG(); 4641 friction->rho_ice=matpar->GetRhoIce(); 4642 friction->rho_water=matpar->GetRhoWater(); 4643 friction->K=&dragcoefficient_list[0]; 4644 friction->bed=&bed_list[0]; 4645 friction->thickness=&thickness_list[0]; 4646 friction->vx=&vx_list[0]; 4647 friction->vy=&vy_list[0]; 4648 friction->p=drag_p; 4649 friction->q=drag_q; 4650 4651 4652 if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!"); 4653 4654 /*Compute alpha complement: */ 4655 FrictionGetAlphaComplement(&alpha_complement_list[0],friction); 4656 4657 /*Erase friction object: */ 4658 DeleteFriction(&friction); 4659 } 4660 else{ 4661 alpha_complement_list[0]=0; 4662 alpha_complement_list[1]=0; 4663 alpha_complement_list[2]=0; 4660 4664 } 4661 4665 … … 4671 4675 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4672 4676 4673 /*Build alpha_complement_list: */4674 if (drag_type==2){4675 4676 /*Allocate friction object: */4677 Friction* friction=NewFriction();4678 4679 /*Initialize all fields: */4680 friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));4681 strcpy(friction->element_type,"2d");4682 4683 friction->gravity=matpar->GetG();4684 friction->rho_ice=matpar->GetRhoIce();4685 friction->rho_water=matpar->GetRhoWater();4686 friction->K=&this->properties.k[0];4687 friction->bed=&this->properties.b[0];4688 friction->thickness=&this->properties.h[0];4689 friction->velocities=&vxvyvz_list[0][0];4690 friction->p=this->properties.p;4691 friction->q=this->properties.q;4692 4693 if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");4694 4695 /*Compute alpha complement: */4696 FrictionGetAlphaComplement(&alpha_complement_list[0],friction);4697 4698 /*Erase friction object: */4699 DeleteFriction(&friction);4700 }4701 else{4702 alpha_complement_list[0]=0;4703 alpha_complement_list[1]=0;4704 alpha_complement_list[2]=0;4705 }4706 4707 4677 /*Recover alpha_complement and k: */ 4708 4678 GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3); 4709 GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3);4679 inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); 4710 4680 4711 4681 /*recover lambda mu and xi: */ 4712 GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);4713 GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);4714 GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3);4682 inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); 4683 inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); 4684 inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); 4715 4685 4716 4686 /*recover vx vy and vz: */ 4717 GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);4718 GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);4719 GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3);4687 inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); 4688 inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); 4689 inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); 4720 4690 4721 4691 /*Get normal vecyor to the bed */ … … 4736 4706 4737 4707 /*Get k derivative: dk/dx */ 4738 GetParameterDerivativeValue(&dk[0], &this->properties.k[0],&xyz_list[0][0], gauss_l1l2l3);4708 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); 4739 4709 4740 4710 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ … … 4748 4718 4749 4719 //Add regularization term 4750 grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);4720 grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); 4751 4721 4752 4722 //min dampening 4753 if(drag< numpar->cm_mindmp_value){4754 grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];4723 if(drag<cm_mindmp_value){ 4724 grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; 4755 4725 } 4756 4726 4757 4727 //max dampening 4758 if(drag> numpar->cm_maxdmp_value){4759 grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];4728 if(drag>cm_maxdmp_value){ 4729 grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; 4760 4730 } 4761 4731 } … … 4839 4809 4840 4810 /*get thickness and velocity at two segment extremities: */ 4841 GetParameterValue(&h1, &this->properties.h[0],gauss_1); 4842 GetParameterValue(&h2, &this->properties.h[0],gauss_2); 4843 GetParameterValue(&vx1, &vx_list[0],gauss_1); 4844 GetParameterValue(&vy1, &vy_list[0],gauss_1); 4845 GetParameterValue(&vx2, &vx_list[0],gauss_2); 4846 GetParameterValue(&vy2, &vy_list[0],gauss_2); 4847 4811 inputs->GetParameterValue(&h1, &gauss_1[0],ThicknessEnum); 4812 inputs->GetParameterValue(&h2, &gauss_2[0],ThicknessEnum); 4813 inputs->GetParameterValue(&vx1, &gauss_1[0],VxEnum); 4814 inputs->GetParameterValue(&vx2, &gauss_2[0],VxEnum); 4815 inputs->GetParameterValue(&vy1, &gauss_1[0],VyEnum); 4816 inputs->GetParameterValue(&vy2, &gauss_2[0],VyEnum); 4848 4817 4849 4818 mass_flux= rho_ice*length*( … … 4871 4840 4872 4841 /* grid data: */ 4873 double vxvy_list[numgrids][2];4874 4842 double vx_list[numgrids]; 4875 4843 double vy_list[numgrids]; 4876 double obs_vxvy_list[numgrids][2];4877 4844 double obs_vx_list[numgrids]; 4878 4845 double obs_vy_list[numgrids]; … … 4888 4855 double gauss_weight; 4889 4856 double gauss_l1l2l3[3]; 4857 double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; 4890 4858 4891 4859 /* parameters: */ … … 4906 4874 Matpar* matpar=NULL; 4907 4875 Matice* matice=NULL; 4908 Numpar* numpar=NULL;4909 4876 4910 4877 /*inputs: */ 4911 4878 bool onwater; 4879 4880 double meanvel, epsvel; 4912 4881 4913 4882 /*retrieve inputs :*/ … … 4921 4890 matpar=(Matpar*)hmatpar.delivers(); 4922 4891 matice=(Matice*)hmatice.delivers(); 4923 numpar=(Numpar*)hnumpar.delivers();4924 4892 4925 4893 /* Get node coordinates and dof list: */ … … 4927 4895 4928 4896 /* Recover input data: */ 4929 if(!inputs->Recover("fit",&fit)) ISSMERROR(" missing fit input parameter"); 4930 if(fit==3 && !inputs->Recover("surfacearea",&S)){ 4931 ISSMERROR("missing surface area input parameter"); 4932 } 4933 if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ 4934 ISSMERROR("missing velocity_obs input parameter"); 4935 } 4936 if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ 4937 ISSMERROR("missing velocity input parameter"); 4938 } 4939 if(!inputs->Recover("weights",&weights_list[0],1,dofs1,numgrids,(void**)nodes)){ 4940 ISSMERROR("missing weights input parameter"); 4941 } 4942 4943 /*Initialize velocities: */ 4944 for(i=0;i<numgrids;i++){ 4945 obs_vx_list[i]=obs_vxvy_list[i][0]; 4946 obs_vy_list[i]=obs_vxvy_list[i][1]; 4947 vx_list[i]=vxvy_list[i][0]; 4948 vy_list[i]=vxvy_list[i][1]; 4949 } 4950 4897 inputs->GetParameterValue(&fit,FitEnum); 4898 if(fit==3){ 4899 inputs->GetParameterValue(&S,SurfaceAreaEnum); 4900 } 4901 inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); 4902 inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); 4903 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); 4904 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); 4905 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 4906 4907 /*retrieve some parameters: */ 4908 this->parameters->FindParam(&meanvel,"meanvel"); 4909 this->parameters->FindParam(&epsvel,"epsvel"); 4910 4951 4911 /* Compute Misfit at the 3 nodes 4952 4912 * Here we integrate linearized functions: … … 4979 4939 */ 4980 4940 for (i=0;i<numgrids;i++){ 4981 scalex=pow( numpar->meanvel/(obs_vx_list[i]+numpar->epsvel),(double)2);4982 scaley=pow( numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),(double)2);4941 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2); 4942 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2); 4983 4943 if(obs_vx_list[i]==0)scalex=0; 4984 4944 if(obs_vy_list[i]==0)scaley=0; … … 4995 4955 */ 4996 4956 for (i=0;i<numgrids;i++){ 4997 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+ numpar->epsvel; //epsvel to avoid velocity being nil.4998 obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+ numpar->epsvel; //epsvel to avoid observed velocity being nil.4999 misfit_list[i]=4*pow( numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2);4957 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil. 4958 obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil. 4959 misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); 5000 4960 } 5001 4961 } … … 5020 4980 */ 5021 4981 for (i=0;i<numgrids;i++){ 5022 misfit_list[i]=0.5*pow( numpar->meanvel,(double)2)*(5023 pow(log((fabs(vx_list[i])+ numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)),(double)2) +5024 pow(log((fabs(vy_list[i])+ numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)),(double)2) );4982 misfit_list[i]=0.5*pow(meanvel,(double)2)*( 4983 pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) + 4984 pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) ); 5025 4985 } 5026 4986 } -
issm/trunk/src/c/objects/Tria.h
r3612 r3620 31 31 Hook hmatice; //hook to 1 matice 32 32 Hook hmatpar; //hook to 1 matpar 33 Hook hnumpar; //hook to 1 numpar34 33 35 Inputs* inputs; 34 Parameters* parameters; //pointer to solution parameters 35 Inputs* inputs; 36 36 37 37 public: … … 39 39 /*FUNCTION constructors, destructors {{{1*/ 40 40 Tria(); 41 Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id , int tria_numpar_id);42 Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs);41 Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id); 42 Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* parameters, Inputs* tria_inputs); 43 43 Tria(int i, IoModel* iomodel); 44 44 ~Tria(); 45 45 /*}}}*/ 46 46 /*FUNCTION object management {{{1*/ 47 void Configure( void* loads,void* nodes,void* materials,void* parameters);47 void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); 48 48 Object* copy(); 49 49 void DeepEcho();
Note:
See TracChangeset
for help on using the changeset viewer.