Changeset 13129
- Timestamp:
- 08/22/12 15:38:04 (13 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h
r13119 r13129 282 282 MacAyeal3dIceFrontEnum, 283 283 MaticeEnum, 284 MatdamageiceEnum, 284 285 MatparEnum, 285 286 NodeEnum, -
issm/trunk-jpl/src/c/Makefile.am
r13086 r13129 109 109 ./classes/objects/Materials/Matice.h\ 110 110 ./classes/objects/Materials/Matice.cpp\ 111 ./classes/objects/Materials/Matdamageice.h\ 112 ./classes/objects/Materials/Matdamageice.cpp\ 111 113 ./classes/objects/Materials/Matpar.h\ 112 114 ./classes/objects/Materials/Matpar.cpp\ -
issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp
r13119 r13129 30 30 31 31 this->nodes=NULL; 32 this->mat ice=NULL;32 this->material=NULL; 33 33 this->matpar=NULL; 34 34 this->verticalneighbors=NULL; … … 49 49 Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels) 50 50 :PentaRef(nummodels) 51 ,PentaHook(nummodels,index+1,iomodel) //index+1: mat iceid, iomodel->numberofelements+1: matpar id51 ,PentaHook(nummodels,index+1,iomodel) //index+1: material id, iomodel->numberofelements+1: matpar id 52 52 { //i is the element index 53 53 … … 85 85 /*initialize pointers:*/ 86 86 this->nodes=NULL; 87 this->mat ice=NULL;87 this->material=NULL; 88 88 this->matpar=NULL; 89 89 this->verticalneighbors=NULL; … … 107 107 penta->hnodes=new Hook*[penta->numanalyses]; 108 108 for(i=0;i<penta->numanalyses;i++)penta->hnodes[i]=(Hook*)this->hnodes[i]->copy(); 109 penta->hmat ice=(Hook*)this->hmatice->copy();109 penta->hmaterial=(Hook*)this->hmaterial->copy(); 110 110 penta->hmatpar=(Hook*)this->hmatpar->copy(); 111 111 penta->hneighbors=(Hook*)this->hneighbors->copy(); … … 132 132 penta->nodes=xNew<Node*>(6); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. 133 133 for(i=0;i<6;i++)penta->nodes[i]=this->nodes[i]; 134 penta->mat ice=(Matice*)penta->hmatice->delivers();134 penta->material=(Material*)penta->hmaterial->delivers(); 135 135 penta->matpar=(Matpar*)penta->hmatpar->delivers(); 136 136 penta->verticalneighbors=(Penta**)penta->hneighbors->deliverp(); … … 286 286 /*Compute strain rate viscosity and pressure: */ 287 287 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 288 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);288 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 289 289 pressure_input->GetInputValue(&pressure,gauss); 290 290 … … 353 353 /*Compute strain rate viscosity and pressure: */ 354 354 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 355 mat ice->GetViscosity3d(&viscosity,&epsilon[0]);355 material->GetViscosity3d(&viscosity,&epsilon[0]); 356 356 pressure_input->GetInputValue(&pressure,gauss); 357 357 … … 391 391 * datasets, using internal ids and offsets hidden in hooks: */ 392 392 if (this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); 393 this->hmat ice->configure(materialsin);393 this->hmaterial->configure(materialsin); 394 394 this->hmatpar->configure(materialsin); 395 395 this->hneighbors->configure(elementsin); … … 398 398 if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 399 399 else this->nodes=NULL; 400 this->mat ice=(Matice*)this->hmatice->delivers();400 this->material=(Material*)this->hmaterial->delivers(); 401 401 this->matpar=(Matpar*)this->hmatpar->delivers(); 402 402 this->verticalneighbors=(Penta**)this->hneighbors->deliverp(); … … 419 419 420 420 /*Checks in debugging {{{*/ 421 _assert_(this->nodes && this->mat ice&& this->matpar && this->verticalneighbors && this->parameters && this->inputs);421 _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); 422 422 /*}}}*/ 423 423 … … 490 490 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 491 491 ElementMatrix* Ke=tria->CreateKMatrixPrognostic(); 492 delete tria->mat ice; delete tria;492 delete tria->material; delete tria; 493 493 494 494 /*Delete Vx and Vy averaged*/ … … 507 507 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 508 508 ElementMatrix* Ke=tria->CreateMassMatrix(); 509 delete tria->mat ice; delete tria;509 delete tria->material; delete tria; 510 510 511 511 /*clean up and return*/ … … 522 522 523 523 /*if debugging mode, check that all pointers exist {{{*/ 524 _assert_(this->nodes && this->mat ice&& this->matpar && this->verticalneighbors && this->parameters && this->inputs);524 _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); 525 525 /*}}}*/ 526 526 … … 591 591 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 592 592 ElementVector* pe=tria->CreatePVectorPrognostic(); 593 delete tria->mat ice; delete tria;593 delete tria->material; delete tria; 594 594 595 595 /*Delete Vx and Vy averaged*/ … … 609 609 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 610 610 ElementVector* pe=tria->CreatePVectorSlope(); 611 delete tria->mat ice; delete tria;611 delete tria->material; delete tria; 612 612 613 613 /*clean up and return*/ … … 624 624 625 625 /*Checks in debugging {{{*/ 626 _assert_(this->nodes && this->mat ice&& this->matpar && this->verticalneighbors && this->parameters && this->inputs);626 _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); 627 627 /*}}}*/ 628 628 … … 661 661 nodes[4]->DeepEcho(); 662 662 nodes[5]->DeepEcho(); 663 mat ice->DeepEcho();663 material->DeepEcho(); 664 664 matpar->DeepEcho(); 665 665 _printLine_(" neighbor ids: " << verticalneighbors[0]->Id() << "-" << verticalneighbors[1]->Id()); … … 1388 1388 original_input=(Input*)penta->inputs->GetInput(enum_type); 1389 1389 else if (object_enum==MaterialsEnum) 1390 original_input=(Input*)penta->mat ice->inputs->GetInput(enum_type);1390 original_input=(Input*)penta->material->inputs->GetInput(enum_type); 1391 1391 else 1392 1392 _error_("object " << EnumToStringx(object_enum) << " not supported yet"); … … 1448 1448 this->inputs->AddInput((Input*)depth_averaged_input); 1449 1449 else if (object_enum==MaterialsEnum) 1450 this->mat ice->inputs->AddInput((Input*)depth_averaged_input);1450 this->material->inputs->AddInput((Input*)depth_averaged_input); 1451 1451 else 1452 1452 _error_("object " << EnumToStringx(object_enum) << " not supported yet"); … … 1481 1481 num_inputs=1; 1482 1482 base_inputs=xNew<Input*>(num_inputs); 1483 base_inputs[0]=(Input*)mat ice->inputs->GetInput(enum_type);1483 base_inputs[0]=(Input*)material->inputs->GetInput(enum_type); 1484 1484 } 1485 1485 else if (object_type==NodeEnum){ … … 1515 1515 } 1516 1516 else if(object_type==MaterialsEnum){ 1517 penta->mat ice->inputs->AddInput((Input*)copy);1517 penta->material->inputs->AddInput((Input*)copy); 1518 1518 } 1519 1519 else if(object_type==NodeEnum){ … … 1554 1554 1555 1555 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ 1556 if (enum_type==MaterialsRheologyBbarEnum) input=this->mat ice->inputs->GetInput(MaterialsRheologyBEnum);1557 else if (enum_type==MaterialsRheologyZbarEnum) input=this->mat ice->inputs->GetInput(MaterialsRheologyZEnum);1556 if (enum_type==MaterialsRheologyBbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyBEnum); 1557 else if (enum_type==MaterialsRheologyZbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyZEnum); 1558 1558 else input=this->inputs->GetInput(enum_type); 1559 1559 //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in penta->inputs"); why error out? if the requested input does not exist, we should still … … 1669 1669 case MaterialsRheologyBbarEnum: 1670 1670 case MaterialsRheologyZbarEnum: 1671 /*Mat icewill take care of it*/ break;1671 /*Material will take care of it*/ break; 1672 1672 default: 1673 1673 _error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); … … 1965 1965 /*update input*/ 1966 1966 if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ 1967 mat ice->inputs->AddInput(new PentaP1Input(name,values));1967 material->inputs->AddInput(new PentaP1Input(name,values)); 1968 1968 } 1969 1969 else{ … … 2637 2637 this->SpawnTriaHook(dynamic_cast<TriaHook*>(tria),&indices[0]); 2638 2638 2639 /*Spawn mat ice*/2640 tria->mat ice=NULL;2641 tria->mat ice=(Matice*)this->matice->copy();2642 delete tria->mat ice->inputs;2643 tria->mat ice->inputs=(Inputs*)this->matice->inputs->SpawnTriaInputs(indices);2644 2645 /*recover nodes, mat iceand matpar: */2639 /*Spawn material*/ 2640 tria->material=NULL; 2641 tria->material=(Material*)this->material->copy(); 2642 delete tria->material->inputs; 2643 tria->material->inputs=(Inputs*)this->material->inputs->SpawnTriaInputs(indices); 2644 2645 /*recover nodes, material and matpar: */ 2646 2646 tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp(); 2647 2647 tria->matpar=(Matpar*)tria->hmatpar->delivers(); … … 2730 2730 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 2731 2731 S=tria->SurfaceArea(); 2732 delete tria->mat ice; delete tria;2732 delete tria->material; delete tria; 2733 2733 return S; 2734 2734 } … … 2737 2737 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 2738 2738 S=tria->SurfaceArea(); 2739 delete tria->mat ice; delete tria;2739 delete tria->material; delete tria; 2740 2740 return S; 2741 2741 } … … 3018 3018 3019 3019 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 3020 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);3020 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 3021 3021 GetPhi(&phi, &epsilon[0], viscosity); 3022 3022 … … 3128 3128 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 3129 3129 mass_flux=tria->MassFlux(segment,process_units); 3130 delete tria->mat ice; delete tria;3130 delete tria->material; delete tria; 3131 3131 3132 3132 /*Delete Vx and Vy averaged*/ … … 3235 3235 switch(response_enum){ 3236 3236 case MaterialsRheologyBbarEnum: 3237 *presponse=this->mat ice->GetBbar();3237 *presponse=this->material->GetBbar(); 3238 3238 break; 3239 3239 case MaterialsRheologyZbarEnum: 3240 *presponse=this->mat ice->GetZbar();3240 *presponse=this->material->GetZbar(); 3241 3241 break; 3242 3242 case VelEnum: … … 3531 3531 ElementMatrix* Ke=tria->CreateKMatrixMelting(); 3532 3532 3533 delete tria->mat ice; delete tria;3533 delete tria->material; delete tria; 3534 3534 return Ke; 3535 3535 } … … 3830 3830 3831 3831 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 3832 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);3832 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 3833 3833 GetPhi(&phi, &epsilon[0], viscosity); 3834 3834 … … 4086 4086 4087 4087 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 4088 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);4088 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 4089 4089 GetPhi(&phi, &epsilon[0], viscosity); 4090 4090 … … 4356 4356 B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); 4357 4357 for(i=0;i<numdof;i++) B[i]=B_average; 4358 this->mat ice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));4358 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); 4359 4359 break; 4360 4360 case ArrheniusEnum: … … 4362 4362 B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, 4363 4363 s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0), 4364 mat ice->GetN());4364 material->GetN()); 4365 4365 for(i=0;i<numdof;i++) B[i]=B_average; 4366 this->mat ice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));4366 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); 4367 4367 break; 4368 4368 default: … … 4434 4434 B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0); 4435 4435 for(i=0;i<numdof;i++) B[i]=B_average; 4436 this->mat ice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));4436 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); 4437 4437 break; 4438 4438 case ArrheniusEnum: … … 4440 4440 B_average=Arrhenius((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0, 4441 4441 s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0), 4442 mat ice->GetN());4442 material->GetN()); 4443 4443 for(i=0;i<numdof;i++) B[i]=B_average; 4444 this->mat ice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));4444 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); 4445 4445 break; 4446 4446 default: … … 4468 4468 if(enum_type==MaterialsRheologyBbarEnum){ 4469 4469 if(!IsOnBed()) return; 4470 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyBEnum);4470 input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); 4471 4471 } 4472 4472 else if(enum_type==MaterialsRheologyZbarEnum){ 4473 4473 if(!IsOnBed()) return; 4474 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyZEnum);4474 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); 4475 4475 } 4476 4476 … … 4491 4491 4492 4492 if(enum_type==MaterialsRheologyBbarEnum){ 4493 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyBEnum);4493 input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); 4494 4494 } 4495 4495 else if(enum_type==MaterialsRheologyZbarEnum){ 4496 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyZEnum);4496 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); 4497 4497 } 4498 4498 else{ … … 4513 4513 4514 4514 if(enum_type==MaterialsRheologyBbarEnum){ 4515 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyBEnum);4515 input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); 4516 4516 } 4517 4517 else if(enum_type==MaterialsRheologyZbarEnum){ 4518 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyZEnum);4518 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); 4519 4519 } 4520 4520 else{ … … 4558 4558 if (!IsOnBed()) return NULL; 4559 4559 4560 /*Depth Averaging B*/ 4561 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 4562 4563 /*Depth Averaging Z*/ 4564 this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); 4560 /*Depth average some fields*/ 4561 switch(this->material->ObjectEnum()){ 4562 case MaticeEnum: 4563 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 4564 break; 4565 case MatdamageiceEnum: 4566 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 4567 this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); 4568 break; 4569 default: 4570 _error_("material "<<EnumToStringx(this->material->ObjectEnum())<<" not supported"); 4571 } 4565 4572 4566 4573 /*Call Tria function*/ 4567 4574 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 4568 4575 ElementMatrix* Ke=tria->CreateKMatrixAdjointMacAyeal(); 4569 delete tria->matice; delete tria; 4570 4571 /*Delete B averaged*/ 4572 this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); 4573 4574 /*Delete Z averaged*/ 4575 this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); 4576 delete tria->material; delete tria; 4577 4578 /*Delete averaged fields*/ 4579 this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); 4580 this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); 4576 4581 4577 4582 /*clean up and return*/ … … 4619 4624 4620 4625 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 4621 mat ice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);4626 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); 4622 4627 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; 4623 4628 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; … … 4689 4694 4690 4695 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 4691 mat ice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);4696 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); 4692 4697 eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; 4693 4698 eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; … … 4754 4759 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 4755 4760 ElementVector* pe=tria->CreatePVectorAdjointHoriz(); 4756 delete tria->mat ice; delete tria;4761 delete tria->material; delete tria; 4757 4762 4758 4763 /*clean up and return*/ … … 4768 4773 Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 4769 4774 ElementVector* pe=tria->CreatePVectorAdjointHoriz(); 4770 delete tria->mat ice; delete tria;4775 delete tria->material; delete tria; 4771 4776 4772 4777 /*clean up and return*/ … … 4782 4787 Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 4783 4788 ElementVector* pe=tria->CreatePVectorAdjointStokes(); 4784 delete tria->mat ice; delete tria;4789 delete tria->material; delete tria; 4785 4790 4786 4791 /*clean up and return*/ … … 4879 4884 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 4880 4885 tria->GradjDragGradient(gradient,resp,control_index); 4881 delete tria->mat ice; delete tria;4886 delete tria->material; delete tria; 4882 4887 break; 4883 4888 case RheologyBbarAbsGradientEnum: … … 4885 4890 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 4886 4891 tria->GradjBGradient(gradient,resp,control_index); 4887 delete tria->mat ice; delete tria;4892 delete tria->material; delete tria; 4888 4893 break; 4889 4894 default: … … 4902 4907 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 4903 4908 tria->GradjDragMacAyeal(gradient,control_index); 4904 delete tria->mat ice; delete tria;4909 delete tria->material; delete tria; 4905 4910 4906 4911 } /*}}}*/ … … 5080 5085 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 5081 5086 tria->GradjBMacAyeal(gradient,control_index); 5082 delete tria->mat ice; delete tria;5087 delete tria->material; delete tria; 5083 5088 5084 5089 /*delete Average B*/ 5085 this->mat ice->inputs->DeleteInput(MaterialsRheologyBbarEnum);5090 this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); 5086 5091 5087 5092 } /*}}}*/ … … 5098 5103 Tria* tria=(Tria*)SpawnTria(0,1,2); 5099 5104 tria->GradjBMacAyeal(gradient,control_index); //We use MacAyeal as an estimate for now 5100 delete tria->mat ice; delete tria;5105 delete tria->material; delete tria; 5101 5106 5102 5107 /*delete Average B*/ 5103 this->mat ice->inputs->DeleteInput(MaterialsRheologyBbarEnum);5108 this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); 5104 5109 } /*}}}*/ 5105 5110 /*FUNCTION Penta::GradjBbarStokes {{{*/ … … 5115 5120 Tria* tria=(Tria*)SpawnTria(0,1,2); 5116 5121 tria->GradjBMacAyeal(gradient,control_index); //We use MacAyeal as an estimate for now 5117 delete tria->mat ice; delete tria;5122 delete tria->material; delete tria; 5118 5123 5119 5124 /*delete Average B*/ 5120 this->mat ice->inputs->DeleteInput(MaterialsRheologyBbarEnum);5125 this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); 5121 5126 } /*}}}*/ 5122 5127 /*FUNCTION Penta::InputControlUpdate{{{*/ … … 5136 5141 if(control_type[i]==MaterialsRheologyBbarEnum){ 5137 5142 if (!IsOnBed()) goto cleanup_and_return; 5138 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);5143 input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input); 5139 5144 } 5140 5145 else if(control_type[i]==MaterialsRheologyZbarEnum){ 5141 5146 if (!IsOnBed()) goto cleanup_and_return; 5142 input=(Input*)mat ice->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input);5147 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input); 5143 5148 } 5144 5149 else{ … … 5268 5273 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 5269 5274 J=tria->SurfaceAverageVelMisfit(process_units,weight_index); 5270 delete tria->mat ice; delete tria;5275 delete tria->material; delete tria; 5271 5276 return J; 5272 5277 } … … 5275 5280 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 5276 5281 J=tria->SurfaceAverageVelMisfit(process_units,weight_index); 5277 delete tria->mat ice; delete tria;5282 delete tria->material; delete tria; 5278 5283 return J; 5279 5284 } … … 5305 5310 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 5306 5311 J=tria->SurfaceAbsVelMisfit(process_units,weight_index); 5307 delete tria->mat ice; delete tria;5312 delete tria->material; delete tria; 5308 5313 return J; 5309 5314 } … … 5312 5317 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 5313 5318 J=tria->SurfaceAbsVelMisfit(process_units,weight_index); 5314 delete tria->mat ice; delete tria;5319 delete tria->material; delete tria; 5315 5320 return J; 5316 5321 } … … 5342 5347 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 5343 5348 J=tria->SurfaceLogVelMisfit(process_units,weight_index); 5344 delete tria->mat ice; delete tria;5349 delete tria->material; delete tria; 5345 5350 return J; 5346 5351 } … … 5349 5354 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 5350 5355 J=tria->SurfaceLogVelMisfit(process_units,weight_index); 5351 delete tria->mat ice; delete tria;5356 delete tria->material; delete tria; 5352 5357 return J; 5353 5358 } … … 5381 5386 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 5382 5387 J=tria->SurfaceLogVxVyMisfit(process_units,weight_index); 5383 delete tria->mat ice; delete tria;5388 delete tria->material; delete tria; 5384 5389 return J; 5385 5390 } … … 5388 5393 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 5389 5394 J=tria->SurfaceLogVxVyMisfit(process_units,weight_index); 5390 delete tria->mat ice; delete tria;5395 delete tria->material; delete tria; 5391 5396 return J; 5392 5397 } … … 5418 5423 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 5419 5424 J=tria->SurfaceRelVelMisfit(process_units,weight_index); 5420 delete tria->mat ice; delete tria;5425 delete tria->material; delete tria; 5421 5426 return J; 5422 5427 } … … 5425 5430 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 5426 5431 J=tria->SurfaceRelVelMisfit(process_units,weight_index); 5427 delete tria->mat ice; delete tria;5432 delete tria->material; delete tria; 5428 5433 return J; 5429 5434 } … … 5452 5457 tria=(Tria*)SpawnTria(0,1,2); 5453 5458 J=tria->ThicknessAbsMisfit(process_units,weight_index); 5454 delete tria->mat ice; delete tria;5459 delete tria->material; delete tria; 5455 5460 return J; 5456 5461 } … … 5467 5472 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria 5468 5473 J=tria->DragCoefficientAbsGradient(process_units,weight_index); 5469 delete tria->mat ice; delete tria;5474 delete tria->material; delete tria; 5470 5475 return J; 5471 5476 } … … 5482 5487 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria 5483 5488 J=tria->RheologyBbarAbsGradient(process_units,weight_index); 5484 delete tria->mat ice; delete tria;5489 delete tria->material; delete tria; 5485 5490 return J; 5486 5491 } … … 5531 5536 5532 5537 if(control_enum==MaterialsRheologyBbarEnum){ 5533 input=(Input*)mat ice->inputs->GetInput(control_enum); _assert_(input);5538 input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); 5534 5539 } 5535 5540 else{ … … 5826 5831 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 5827 5832 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 5828 mat ice->GetViscosity3d(&viscosity, &epsilon[0]);5829 mat ice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);5833 material->GetViscosity3d(&viscosity, &epsilon[0]); 5834 material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); 5830 5835 5831 5836 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); … … 5847 5852 5848 5853 /*Clean-up and return*/ 5849 delete tria->mat ice; delete tria;5854 delete tria->material; delete tria; 5850 5855 delete gauss; 5851 5856 delete gauss_tria; … … 6036 6041 6037 6042 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 6038 mat ice->GetViscosity3dStokes(&viscosity, &epsilon[0]);6043 material->GetViscosity3dStokes(&viscosity, &epsilon[0]); 6039 6044 6040 6045 D_scalar=2*viscosity*gauss->weight*Jdet; … … 6063 6068 6064 6069 /*Clean-up and return*/ 6065 delete tria->mat ice; delete tria;6070 delete tria->material; delete tria; 6066 6071 delete gauss; 6067 6072 delete gauss_tria; … … 6143 6148 6144 6149 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 6145 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);6150 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 6146 6151 6147 6152 BedNormal(&bed_normal[0],xyz_list_tria); … … 6335 6340 if (!IsOnBed()) return NULL; 6336 6341 6337 /*Depth Averaging B*/ 6338 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 6339 6340 /*Depth Averaging Z*/ 6341 this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); 6342 /*Depth average some fields*/ 6343 switch(this->material->ObjectEnum()){ 6344 case MaticeEnum: 6345 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 6346 break; 6347 case MatdamageiceEnum: 6348 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 6349 this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); 6350 break; 6351 default: 6352 _error_("material "<<EnumToStringx(this->material->ObjectEnum())<<" not supported"); 6353 } 6342 6354 6343 6355 /*Call Tria function*/ 6344 6356 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 6345 6357 ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyeal(); 6346 delete tria->matice; delete tria; 6347 6348 /*Delete B averaged*/ 6349 this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); 6350 6351 /*Delete Z averaged*/ 6352 this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); 6358 delete tria->material; delete tria; 6359 6360 /*Delete averaged fields*/ 6361 this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); 6362 this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); 6353 6363 6354 6364 /*clean up and return*/ … … 6425 6435 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 6426 6436 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 6427 mat ice->GetViscosity3d(&viscosity, &epsilon[0]);6428 mat ice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);6437 material->GetViscosity3d(&viscosity, &epsilon[0]); 6438 material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); 6429 6439 6430 6440 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); … … 6432 6442 else if (approximation==MacAyealStokesApproximationEnum){ 6433 6443 this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 6434 mat ice->GetViscosity3dStokes(&newviscosity,&epsilons[0]);6444 material->GetViscosity3dStokes(&newviscosity,&epsilons[0]); 6435 6445 } 6436 6446 else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); … … 6451 6461 6452 6462 /*Clean up and return*/ 6453 delete tria->mat ice;6463 delete tria->material; 6454 6464 delete tria; 6455 6465 delete gauss_tria; … … 6469 6479 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 6470 6480 ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction(); 6471 delete tria->mat ice; delete tria;6481 delete tria->material; delete tria; 6472 6482 6473 6483 /*clean-up and return*/ … … 6582 6592 6583 6593 /*Clean up and return*/ 6584 delete tria->mat ice;6594 delete tria->material; 6585 6595 delete tria; 6586 6596 delete gauss_tria; … … 6600 6610 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 6601 6611 ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction(); 6602 delete tria->mat ice; delete tria;6612 delete tria->material; delete tria; 6603 6613 6604 6614 /*clean-up and return*/ … … 6665 6675 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 6666 6676 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 6667 mat ice->GetViscosity3d(&viscosity, &epsilon[0]);6668 mat ice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);6677 material->GetViscosity3d(&viscosity, &epsilon[0]); 6678 material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); 6669 6679 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 6670 6680 … … 6829 6839 6830 6840 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 6831 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);6841 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 6832 6842 6833 6843 D_scalar=gauss->weight*Jdet; … … 6900 6910 6901 6911 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 6902 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);6912 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 6903 6913 6904 6914 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); … … 7088 7098 7089 7099 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 7090 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);7100 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 7091 7101 7092 7102 for(i=0;i<NUMVERTICES;i++){ … … 7162 7172 BedNormal(&bed_normal[0],xyz_list_tria); 7163 7173 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 7164 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);7174 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 7165 7175 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); 7166 7176 … … 7239 7249 7240 7250 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 7241 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);7251 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 7242 7252 7243 7253 for(i=0;i<NUMVERTICES;i++){ … … 7313 7323 BedNormal(&bed_normal[0],xyz_list_tria); 7314 7324 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 7315 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);7325 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 7316 7326 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); 7317 7327 … … 7435 7445 rho_ice=matpar->GetRhoIce(); 7436 7446 gravity=matpar->GetG(); 7437 n=mat ice->GetN();7438 B=mat ice->GetB();7447 n=material->GetN(); 7448 B=material->GetB(); 7439 7449 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); 7440 7450 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); … … 7504 7514 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 7505 7515 ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); 7506 delete tria->mat ice; delete tria;7516 delete tria->material; delete tria; 7507 7517 7508 7518 /*Clean up and return*/ … … 7518 7528 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 7519 7529 ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); 7520 delete tria->mat ice; delete tria;7530 delete tria->material; delete tria; 7521 7531 7522 7532 /*Clean up and return*/ … … 7636 7646 7637 7647 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 7638 mat ice->GetViscosity3dStokes(&viscosity,&epsilon[0]);7648 material->GetViscosity3dStokes(&viscosity,&epsilon[0]); 7639 7649 7640 7650 for(i=0;i<NUMVERTICES+1;i++){ … … 7897 7907 if (!IsOnBed()) return NULL; 7898 7908 7899 /*Depth Averaging B*/ 7900 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 7901 7902 /*Depth Averaging Z*/ 7903 this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); 7909 /*Depth average some fields*/ 7910 switch(this->material->ObjectEnum()){ 7911 case MaticeEnum: 7912 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 7913 break; 7914 case MatdamageiceEnum: 7915 this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); 7916 this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); 7917 break; 7918 default: 7919 _error_("material "<<EnumToStringx(this->material->ObjectEnum())<<" not supported"); 7920 } 7904 7921 7905 7922 /*Call Tria function*/ 7906 7923 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 7907 7924 ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal(); 7908 delete tria->matice; delete tria; 7909 7910 /*Delete B averaged*/ 7911 this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); 7912 7913 /*Delete Z averaged*/ 7914 this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); 7925 delete tria->material; delete tria; 7926 7927 /*Delete averaged inputs*/ 7928 this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); 7929 this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); 7915 7930 7916 7931 /*clean up and return*/ … … 7955 7970 7956 7971 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 7957 mat ice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);7972 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); 7958 7973 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; 7959 7974 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; … … 8022 8037 8023 8038 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 8024 mat ice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);8039 material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); 8025 8040 eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; 8026 8041 eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; … … 8264 8279 8265 8280 /*Get A*/ 8266 _assert_(mat ice->GetN()==3.0);8267 A=mat ice->GetA();8281 _assert_(material->GetN()==3.0); 8282 A=material->GetA(); 8268 8283 8269 8284 /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/ … … 9085 9100 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 9086 9101 ElementMatrix* Ke=tria->CreateKMatrixBalancethickness(); 9087 delete tria->mat ice; delete tria;9102 delete tria->material; delete tria; 9088 9103 9089 9104 /*Delete Vx and Vy averaged*/ … … 9107 9122 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 9108 9123 ElementVector* pe=tria->CreatePVectorBalancethickness(); 9109 delete tria->mat ice; delete tria;9124 delete tria->material; delete tria; 9110 9125 9111 9126 /*Delete Vx and Vy averaged*/ -
issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h
r13092 r13129 16 16 class IoModel; 17 17 class Node; 18 class Mat ice;18 class Material; 19 19 class Matpar; 20 20 class Tria; … … 35 35 36 36 Node **nodes; // 6 nodes 37 Mat ice *matice; // 1 material ice37 Material *material; // 1 material ice 38 38 Matpar *matpar; // 1 material parameter 39 39 Penta **verticalneighbors; // 2 neighbors: first one under, second one above -
issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.cpp
r12832 r13129 25 25 numanalyses=UNDEF; 26 26 this->hnodes=NULL; 27 this->hmat ice=NULL;27 this->hmaterial=NULL; 28 28 this->hmatpar=NULL; 29 29 this->hneighbors=NULL; … … 39 39 } 40 40 delete [] this->hnodes; 41 delete hmat ice;41 delete hmaterial; 42 42 delete hmatpar; 43 43 delete hneighbors; 44 44 } 45 45 /*}}}*/ 46 /*FUNCTION PentaHook::PentaHook(int in_numanalyses,int mat ice_id, int matpar_id){{{*/47 PentaHook::PentaHook(int in_numanalyses,int mat ice_id, IoModel* iomodel){46 /*FUNCTION PentaHook::PentaHook(int in_numanalyses,int material_id, int matpar_id){{{*/ 47 PentaHook::PentaHook(int in_numanalyses,int material_id, IoModel* iomodel){ 48 48 49 49 /*intermediary: */ … … 55 55 this->numanalyses=in_numanalyses; 56 56 this->hnodes=new Hook*[in_numanalyses]; 57 this->hmat ice=new Hook(&matice_id,1);57 this->hmaterial=new Hook(&material_id,1); 58 58 this->hmatpar=new Hook(&matpar_id,1); 59 59 this->hneighbors=NULL; … … 98 98 } 99 99 } 100 // do not spawn hmat ice. maticewill be taken care of by Penta101 triahook->hmat ice=NULL;100 // do not spawn hmaterial. material will be taken care of by Penta 101 triahook->hmaterial=NULL; 102 102 triahook->hmatpar=(Hook*)this->hmatpar->copy(); 103 103 } -
issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.h
r12365 r13129 15 15 int numanalyses; //number of analysis types 16 16 Hook** hnodes; // 6 nodes for each analysis type 17 Hook* hmat ice; // 1 ice material17 Hook* hmaterial; // 1 ice material 18 18 Hook* hmatpar; // 1 material parameter 19 19 Hook* hneighbors; // 2 elements, first down, second up … … 21 21 /*FUNCTION constructors, destructors {{{*/ 22 22 PentaHook(); 23 PentaHook(int in_numanalyses,int mat ice_id, IoModel* iomodel);23 PentaHook(int in_numanalyses,int material_id, IoModel* iomodel); 24 24 ~PentaHook(); 25 25 void SetHookNodes(int* node_ids,int analysis_counter); -
issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
r13119 r13129 29 29 30 30 this->nodes=NULL; 31 this->mat ice=NULL;31 this->material=NULL; 32 32 this->matpar=NULL; 33 33 for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF; … … 61 61 /*initialize pointers:*/ 62 62 this->nodes=NULL; 63 this->mat ice=NULL;63 this->material=NULL; 64 64 this->matpar=NULL; 65 65 … … 89 89 tria->hnodes=new Hook*[tria->numanalyses]; 90 90 for(i=0;i<tria->numanalyses;i++)tria->hnodes[i]=(Hook*)this->hnodes[i]->copy(); 91 tria->hmat ice=(Hook*)this->hmatice->copy();91 tria->hmaterial=(Hook*)this->hmaterial->copy(); 92 92 tria->hmatpar=(Hook*)this->hmatpar->copy(); 93 93 … … 113 113 tria->nodes=xNew<Node*>(3); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. 114 114 for(i=0;i<3;i++)tria->nodes[i]=this->nodes[i]; 115 tria->mat ice=(Matice*)tria->hmatice->delivers();115 tria->material=(Material*)tria->hmaterial->delivers(); 116 116 tria->matpar=(Matpar*)tria->hmatpar->delivers(); 117 117 … … 173 173 174 174 /*Checks in debugging mode{{{*/ 175 _assert_(this->nodes && this->mat ice&& this->matpar && this->parameters && this->inputs);175 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); 176 176 /*}}}*/ 177 177 … … 526 526 /*asserts: {{{*/ 527 527 /*if debugging mode, check that all pointers exist*/ 528 _assert_(this->nodes && this->mat ice&& this->matpar && this->parameters && this->inputs);528 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); 529 529 /*}}}*/ 530 530 … … 748 748 749 749 /*Checks in debugging {{{*/ 750 _assert_(this->nodes && this->mat ice&& this->matpar && this->parameters && this->inputs);750 _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); 751 751 /*}}}*/ 752 752 … … 812 812 /*Compute strain rate viscosity and pressure: */ 813 813 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 814 mat ice->GetViscosity2d(&viscosity,&epsilon[0]);814 material->GetViscosity2d(&viscosity,&epsilon[0]); 815 815 pressure_input->GetInputValue(&pressure,gauss); 816 816 … … 846 846 * datasets, using internal ids and offsets hidden in hooks: */ 847 847 if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); 848 this->hmat ice->configure(materialsin);848 this->hmaterial->configure(materialsin); 849 849 this->hmatpar->configure(materialsin); 850 850 … … 852 852 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 853 853 else this->nodes=NULL; 854 this->mat ice=(Matice*)this->hmatice->delivers();854 this->material=(Material*)this->hmaterial->delivers(); 855 855 this->matpar=(Matpar*)this->hmatpar->delivers(); 856 856 … … 875 875 else _printLine_("nodes = NULL"); 876 876 877 if (mat ice) matice->DeepEcho();878 else _printLine_("mat ice= NULL");877 if (material) material->DeepEcho(); 878 else _printLine_("material = NULL"); 879 879 880 880 if (matpar) matpar->DeepEcho(); … … 987 987 else _printLine_("nodes = NULL"); 988 988 989 if (mat ice) matice->Echo();990 else _printLine_("mat ice= NULL");989 if (material) material->Echo(); 990 else _printLine_("material = NULL"); 991 991 992 992 if (matpar) matpar->Echo(); … … 1350 1350 oldinput=(Input*)this->inputs->GetInput(enum_type); 1351 1351 else if (object_enum==MaterialsEnum) 1352 oldinput=(Input*)this->mat ice->inputs->GetInput(enum_type);1352 oldinput=(Input*)this->material->inputs->GetInput(enum_type); 1353 1353 else 1354 1354 _error_("object " << EnumToStringx(object_enum) << " not supported yet"); … … 1363 1363 this->inputs->AddInput((Input*)newinput); 1364 1364 else if (object_enum==MaterialsEnum) 1365 this->mat ice->inputs->AddInput((Input*)newinput);1365 this->material->inputs->AddInput((Input*)newinput); 1366 1366 else 1367 1367 _error_("object " << EnumToStringx(object_enum) << " not supported yet"); … … 1396 1396 1397 1397 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ 1398 if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->mat ice->inputs->GetInput(enum_type);1398 if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->material->inputs->GetInput(enum_type); 1399 1399 else input=this->inputs->GetInput(enum_type); 1400 1400 //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs"); … … 1503 1503 case MaterialsRheologyBbarEnum: 1504 1504 case MaterialsRheologyZbarEnum: 1505 /*Mat icewill take care of it*/ break;1505 /*Material will take care of it*/ break; 1506 1506 default: 1507 1507 _error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); … … 1692 1692 /*update input*/ 1693 1693 if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ 1694 mat ice->inputs->AddInput(new TriaP1Input(name,values));1694 material->inputs->AddInput(new TriaP1Input(name,values)); 1695 1695 } 1696 1696 else{ … … 2792 2792 switch(response_enum){ 2793 2793 case MaterialsRheologyBbarEnum: 2794 *presponse=this->mat ice->GetBbar();2794 *presponse=this->material->GetBbar(); 2795 2795 break; 2796 2796 case MaterialsRheologyZbarEnum: 2797 *presponse=this->mat ice->GetZbar();2797 *presponse=this->material->GetZbar(); 2798 2798 break; 2799 2799 case VelEnum:{ … … 2909 2909 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 2910 2910 this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 2911 mat ice->GetViscosity2d(&viscosity, &epsilon[0]);2912 mat ice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);2911 material->GetViscosity2d(&viscosity, &epsilon[0]); 2912 material->GetViscosity2d(&oldviscosity, &oldepsilon[0]); 2913 2913 thickness_input->GetInputValue(&thickness, gauss); 2914 2914 … … 3092 3092 rho_ice=matpar->GetRhoIce(); 3093 3093 gravity=matpar->GetG(); 3094 n=mat ice->GetN();3095 B=mat ice->GetBbar();3094 n=material->GetN(); 3095 B=material->GetBbar(); 3096 3096 Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); 3097 3097 Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); … … 3164 3164 thickness_input->GetInputValue(&thickness, gauss); 3165 3165 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 3166 mat ice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);3166 material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); 3167 3167 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; 3168 3168 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; … … 3409 3409 3410 3410 if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){ 3411 input=(Input*)mat ice->inputs->GetInput(control_type[i]); _assert_(input);3411 input=(Input*)material->inputs->GetInput(control_type[i]); _assert_(input); 3412 3412 } 3413 3413 else{ … … 3436 3436 3437 3437 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ 3438 input=(Input*)mat ice->inputs->GetInput(enum_type);3438 input=(Input*)material->inputs->GetInput(enum_type); 3439 3439 } 3440 3440 else{ … … 3454 3454 3455 3455 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ 3456 input=(Input*)mat ice->inputs->GetInput(enum_type);3456 input=(Input*)material->inputs->GetInput(enum_type); 3457 3457 } 3458 3458 else{ … … 3473 3473 3474 3474 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ 3475 input=(Input*)mat ice->inputs->GetInput(enum_type);3475 input=(Input*)material->inputs->GetInput(enum_type); 3476 3476 } 3477 3477 else{ … … 3567 3567 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 3568 3568 GradientIndexing(&doflist1[0],control_index); 3569 Input* rheologyb_input=mat ice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);3569 Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); 3570 3570 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 3571 3571 … … 3607 3607 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 3608 3608 GradientIndexing(&doflist1[0],control_index); 3609 Input* rheologyz_input=mat ice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);3609 Input* rheologyz_input=material->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); 3610 3610 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 3611 3611 … … 3656 3656 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); 3657 3657 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); 3658 Input* rheologyb_input=mat ice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);3658 Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); 3659 3659 3660 3660 /* Start looping on the number of gaussian points: */ … … 3672 3672 3673 3673 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 3674 mat ice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);3674 material->GetViscosityComplement(&viscosity_complement,&epsilon[0]); 3675 3675 3676 3676 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); … … 3713 3713 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); 3714 3714 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); 3715 Input* rheologyz_input=mat ice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);3715 Input* rheologyz_input=material->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); 3716 3716 3717 3717 /* Start looping on the number of gaussian points: */ … … 3729 3729 3730 3730 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 3731 mat ice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);3731 material->GetViscosityZComplement(&viscosity_complement,&epsilon[0]); 3732 3732 3733 3733 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); … … 4013 4013 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 4014 4014 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 4015 Input* rheologyb_input=mat ice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);4015 Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); 4016 4016 4017 4017 /* Start looping on the number of gaussian points: */ … … 5110 5110 thickness_input->GetInputValue(&thickness, gauss); 5111 5111 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 5112 mat ice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);5112 material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); 5113 5113 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; 5114 5114 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; … … 5216 5216 /*Get input (either in element or material)*/ 5217 5217 if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ 5218 input=(Input*)mat ice->inputs->GetInput(control_enum); _assert_(input);5218 input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); 5219 5219 } 5220 5220 else{ … … 5251 5251 5252 5252 if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ 5253 input=(Input*)mat ice->inputs->GetInput(control_enum); _assert_(input);5253 input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); 5254 5254 } 5255 5255 else{ -
issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h
r13119 r13129 15 15 class IoModel; 16 16 class Node; 17 class Mat ice;17 class Material; 18 18 class Matpar; 19 19 class ElementMatrix; … … 33 33 34 34 Node **nodes; // 3 nodes 35 Mat ice *matice; // 1 material ice35 Material *material; // 1 material ice 36 36 Matpar *matpar; // 1 material parameter 37 37 int horizontalneighborsids[3]; -
issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.cpp
r12832 r13129 25 25 numanalyses=UNDEF; 26 26 this->hnodes=NULL; 27 this->hmat ice=NULL;27 this->hmaterial=NULL; 28 28 this->hmatpar=NULL; 29 29 } … … 37 37 } 38 38 delete [] this->hnodes; 39 delete hmat ice;39 delete hmaterial; 40 40 delete hmatpar; 41 41 42 42 } 43 43 /*}}}*/ 44 /*FUNCTION TriaHook::TriaHook(int in_numanalyses,int mat ice_id, int matpar_id){{{*/45 TriaHook::TriaHook(int in_numanalyses,int mat ice_id, IoModel* iomodel){44 /*FUNCTION TriaHook::TriaHook(int in_numanalyses,int material_id, int matpar_id){{{*/ 45 TriaHook::TriaHook(int in_numanalyses,int material_id, IoModel* iomodel){ 46 46 47 47 /*intermediary: */ … … 53 53 this->numanalyses=in_numanalyses; 54 54 this->hnodes= new Hook*[in_numanalyses]; 55 this->hmat ice=new Hook(&matice_id,1);55 this->hmaterial=new Hook(&material_id,1); 56 56 this->hmatpar=new Hook(&matpar_id,1); 57 57 -
issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.h
r12365 r13129 14 14 int numanalyses; //number of analysis types 15 15 Hook** hnodes; // 3 nodes for each analysis type 16 Hook* hmat ice; // 1 ice material16 Hook* hmaterial; // 1 ice material 17 17 Hook* hmatpar; // 1 material parameter 18 18 … … 20 20 /*FUNCTION constructors, destructors {{{*/ 21 21 TriaHook(); 22 TriaHook(int in_numanalyses,int mat ice_id, IoModel* iomodel);22 TriaHook(int in_numanalyses,int material_id, IoModel* iomodel); 23 23 ~TriaHook(); 24 24 void SetHookNodes(int* node_ids,int analysis_counter); -
issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp
r13036 r13129 498 498 499 499 /*clean-up and return*/ 500 delete tria->mat ice;500 delete tria->material; 501 501 delete tria; 502 502 delete icefront; -
issm/trunk-jpl/src/c/classes/objects/Materials/Material.h
r12832 r13129 17 17 18 18 public: 19 Inputs* inputs; 19 20 virtual ~Material(){}; 21 /*WARNING: input should not be public but it is an easy way to update B from T (using UpdateFromSolution) from Pentas*/ 20 22 21 23 /*Numerics*/ 22 virtual void InputDuplicate(int original_enum,int new_enum)=0; 23 virtual void Configure(Elements* elements)=0; 24 virtual void GetVectorFromInputs(Vector* vector,int input_enum)=0; 24 virtual void InputDuplicate(int original_enum,int new_enum)=0; 25 virtual void Configure(Elements* elements)=0; 26 virtual void GetVectorFromInputs(Vector* vector,int input_enum)=0; 27 virtual void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; 28 virtual void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0; 29 virtual void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0; 30 virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; 31 virtual void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; 32 virtual void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; 33 virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; 34 virtual IssmDouble GetA()=0; 35 virtual IssmDouble GetB()=0; 36 virtual IssmDouble GetBbar()=0; 37 virtual IssmDouble GetN()=0; 38 virtual IssmDouble GetZ()=0; 39 virtual IssmDouble GetZbar()=0; 25 40 26 41 }; -
issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp
r13119 r13129 172 172 } 173 173 /*}}}*/ 174 /*FUNCTION Matice::GetZ {{{*/175 IssmDouble Matice::GetZ(){176 177 /*Output*/178 IssmDouble Z;179 180 inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);181 return Z;182 }183 /*}}}*/184 /*FUNCTION Matice::GetZbar {{{*/185 IssmDouble Matice::GetZbar(){186 187 /*Output*/188 IssmDouble Zbar;189 190 inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);191 return Zbar;192 }193 /*}}}*/194 174 /*FUNCTION Matice::GetVectorFromInputs{{{*/ 195 175 void Matice::GetVectorFromInputs(Vector* vector,int input_enum){ … … 227 207 void Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ 228 208 /*From a string tensor and a material object, return viscosity, using Glen's flow law. 229 Z *B209 B 230 210 viscosity= ------------------------------------------------------------------- 231 211 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] … … 246 226 /*Intermediary: */ 247 227 IssmDouble A,e; 248 IssmDouble B tmp,B,n,Z;228 IssmDouble B,n; 249 229 250 230 /*Get B and n*/ 251 Btmp=GetBbar(); 252 Z=GetZbar(); 231 B=GetBbar(); 253 232 n=GetN(); 254 B=Z*Btmp;255 233 256 234 if (n==1){ … … 314 292 /*Intermediaries: */ 315 293 IssmDouble A,e; 316 IssmDouble B,n,Z; 317 318 /*Get B, Z and n*/ 294 IssmDouble B,n; 295 296 /*Get B and n*/ 297 B=GetB(); 319 298 n=GetN(); 320 Z=GetZ();321 B=Z*GetB();322 299 323 300 if (n==1){ … … 385 362 /*Intermediaries: */ 386 363 IssmDouble A,e; 387 IssmDouble B,n ,Z;364 IssmDouble B,n; 388 365 IssmDouble eps0; 389 366 390 367 /*Get B and n*/ 391 368 eps0=pow((IssmDouble)10,(IssmDouble)-27); 369 B=GetB(); 392 370 n=GetN(); 393 Z=GetZ();394 B=Z*GetB();395 371 396 372 if (n==1){ … … 476 452 477 453 viscosity_complement=1/(2*pow(A,e)); 478 }479 }480 else{481 viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);482 }483 484 /*Checks in debugging mode*/485 _assert_(B>0);486 _assert_(n>0);487 _assert_(viscosity_complement>0);488 489 /*Return: */490 *pviscosity_complement=viscosity_complement;491 }492 /*}}}*/493 /*FUNCTION Matice::GetViscosityZComplement {{{*/494 void Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){495 /*Return viscosity derivative for control method d(mu)/dZ:496 *497 * B498 * dviscosity= -------------------------------------------------------------------499 * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]500 *501 * If epsilon is NULL, it means this is the first time Gradjb is being run, and we502 * return mu20, initial viscosity.503 */504 505 /*output: */506 IssmDouble viscosity_complement;507 508 /*input strain rate: */509 IssmDouble exx,eyy,exy;510 511 /*Intermediary value A and exponent e: */512 IssmDouble A,e;513 IssmDouble B,n;514 515 /*Get B and n*/516 B=GetBbar();517 n=GetN();518 519 if(epsilon){520 exx=*(epsilon+0);521 eyy=*(epsilon+1);522 exy=*(epsilon+2);523 524 /*Build viscosity: mu2=B/(2*A^e) */525 A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;526 if(A==0){527 /*Maximum viscosity_complement for 0 shear areas: */528 viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);529 }530 else{531 e=(n-1)/(2*n);532 533 viscosity_complement=B/(2*pow(A,e));534 454 } 535 455 } … … 781 701 } 782 702 783 /*Get Z*/784 if (iomodel->Data(MaterialsRheologyZEnum)) {785 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];786 this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));787 }788 789 703 /*Control Inputs*/ 790 704 #ifdef _HAVE_CONTROL_ … … 801 715 } 802 716 break; 803 case MaterialsRheologyZbarEnum:804 if (iomodel->Data(MaterialsRheologyZEnum)){805 _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));806 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];807 for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];808 for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];809 this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));810 }811 break;812 813 717 } 814 718 } … … 837 741 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index]; 838 742 this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); 839 }840 841 /*Get Z*/842 if (iomodel->Data(MaterialsRheologyZEnum)) {843 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];844 this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));845 743 } 846 744 … … 859 757 } 860 758 break; 861 case MaterialsRheologyZbarEnum:862 if (iomodel->Data(MaterialsRheologyZEnum)){863 _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));864 for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];865 for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];866 for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];867 this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));868 }869 break;870 759 } 871 760 } … … 886 775 name==MaterialsRheologyBEnum || 887 776 name==MaterialsRheologyBbarEnum || 888 name==MaterialsRheologyNEnum || 889 name==MaterialsRheologyZEnum || 890 name==MaterialsRheologyZbarEnum 777 name==MaterialsRheologyNEnum 891 778 ){ 892 779 return true; -
issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h
r13119 r13129 15 15 16 16 private: 17 /*Id*/18 17 int mid; 19 20 /*hooks: */21 18 Hook* helement; 22 19 23 20 public: 24 /*WARNING: input should not be public but it is an easy way to update B from T (using UpdateFromSolution) from Pentas*/25 Inputs* inputs;26 27 21 /*Matice constructors, destructors: {{{*/ 28 22 Matice(); … … 56 50 void Configure(Elements* elements); 57 51 void GetVectorFromInputs(Vector* vector,int input_enum); 58 /*}}}*/ 59 /*Matice Numerics: {{{*/ 60 void SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); 61 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); 62 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); 63 void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); 64 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); 65 void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); 66 void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 67 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 52 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); 53 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); 54 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); 55 void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); 56 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); 57 void GetViscosityZComplement(IssmDouble*, IssmDouble*){_error_("not supported");}; 58 void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 59 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 68 60 IssmDouble GetA(); 69 61 IssmDouble GetB(); 70 62 IssmDouble GetBbar(); 63 IssmDouble GetZ(){_error_("not supported");}; 64 IssmDouble GetZbar(){_error_("not supported");}; 71 65 IssmDouble GetN(); 72 IssmDouble GetZ(); 73 IssmDouble GetZbar(); 74 bool IsInput(int name); 66 bool IsInput(int name); 75 67 /*}}}*/ 76 68 }; -
issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp
r12832 r13129 45 45 iomodel->Constant(&this->hydro_p,HydrologyPEnum); 46 46 iomodel->Constant(&this->hydro_q,HydrologyQEnum); 47 this->inputs=NULL; /*not used here*/ 47 48 } 48 49 /*}}}*/ -
issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h
r13056 r13129 68 68 void Configure(Elements* elements); 69 69 void GetVectorFromInputs(Vector* vector,int input_enum){return;} 70 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");}; 71 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");}; 72 void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");}; 73 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; 74 void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; 75 void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; 76 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; 77 IssmDouble GetA(){_error_("not supported");}; 78 IssmDouble GetB(){_error_("not supported");}; 79 IssmDouble GetBbar(){_error_("not supported");}; 80 IssmDouble GetN(){_error_("not supported");}; 81 IssmDouble GetZ(){_error_("not supported");}; 82 IssmDouble GetZbar(){_error_("not supported");}; 70 83 /*}}}*/ 71 84 /*Numerics: {{{*/ -
issm/trunk-jpl/src/c/classes/objects/Params/BoolParam.h
r13036 r13129 43 43 int InstanceEnum(){return enum_type;} 44 44 void GetParameterValue(bool* pbool){*pbool=value;} 45 void GetParameterValue(int* pinteger){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return an integer");}46 void GetParameterValue(int** pintarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return an array of integers");}47 void GetParameterValue(int** pintarray,int* pM,int* pN){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return an array of integers");}48 void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a IssmDouble");}45 void GetParameterValue(int* pinteger){_error_("Param "<< EnumToStringx(enum_type) << " cannot return an integer");} 46 void GetParameterValue(int** pintarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return an array of integers");} 47 void GetParameterValue(int** pintarray,int* pM,int* pN){_error_("Param "<< EnumToStringx(enum_type) << " cannot return an array of integers");} 48 void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} 49 49 void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} 50 void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a string");}51 void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a string array");}52 void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a IssmDouble array");}53 void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, int* pN){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a IssmDouble array");}54 void GetParameterValue(IssmDouble*** parray, int* pM,int** pmdims, int** pndims){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a matrix array");}55 void GetParameterValue(Vector** pvec){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a Vec");}56 void GetParameterValue(Matrix** pmat){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a Mat");}57 void GetParameterValue(FILE** pfid){_error_("Param "<< EnumToStringx(enum_type) << " )cannot return a FILE");}50 void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} 51 void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} 52 void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble array");} 53 void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, int* pN){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble array");} 54 void GetParameterValue(IssmDouble*** parray, int* pM,int** pmdims, int** pndims){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a matrix array");} 55 void GetParameterValue(Vector** pvec){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a Vec");} 56 void GetParameterValue(Matrix** pmat){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a Mat");} 57 void GetParameterValue(FILE** pfid){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a FILE");} 58 58 59 59 void SetValue(bool boolean){this->value=boolean;} 60 void SetValue(int integer){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold an int");}61 void SetValue(IssmDouble scalar){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold an IssmPDouble");}62 void SetValue(char* string){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a string");}63 void SetValue(char** stringarray,int M){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a string array");}64 void SetValue(IssmDouble* IssmDoublearray,int M){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a IssmDouble array");}65 void SetValue(IssmDouble* pIssmDoublearray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a IssmDouble array");}66 void SetValue(int* intarray,int M){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a int array");}67 void SetValue(int* pintarray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a int array");}68 void SetValue(Vector* vec){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a Vec");}69 void SetValue(Matrix* mat){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a Mat");}70 void SetValue(FILE* fid){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold a FILE");}71 void SetValue(IssmDouble** array, int M, int* mdim_array, int* ndim_array){_error_("Param "<< EnumToStringx(enum_type) << " )cannot hold an array of matrices");}60 void SetValue(int integer){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an int");} 61 void SetValue(IssmDouble scalar){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an IssmPDouble");} 62 void SetValue(char* string){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a string");} 63 void SetValue(char** stringarray,int M){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a string array");} 64 void SetValue(IssmDouble* IssmDoublearray,int M){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a IssmDouble array");} 65 void SetValue(IssmDouble* pIssmDoublearray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a IssmDouble array");} 66 void SetValue(int* intarray,int M){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a int array");} 67 void SetValue(int* pintarray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a int array");} 68 void SetValue(Vector* vec){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a Vec");} 69 void SetValue(Matrix* mat){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a Mat");} 70 void SetValue(FILE* fid){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a FILE");} 71 void SetValue(IssmDouble** array, int M, int* mdim_array, int* ndim_array){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an array of matrices");} 72 72 void UnitConversion(int direction_enum); 73 73 -
issm/trunk-jpl/src/c/classes/objects/objects.h
r12835 r13129 108 108 #include "./Materials/Material.h" 109 109 #include "./Materials/Matice.h" 110 #include "./Materials/Matdamageice.h" 110 111 #include "./Materials/Matpar.h" 111 112 -
issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp
r13119 r13129 279 279 case MacAyeal3dIceFrontEnum : return "MacAyeal3dIceFront"; 280 280 case MaticeEnum : return "Matice"; 281 case MatdamageiceEnum : return "Matdamageice"; 281 282 case MatparEnum : return "Matpar"; 282 283 case NodeEnum : return "Node"; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
r13119 r13129 17 17 /*Intermediary*/ 18 18 int i,j,k,n; 19 int dim ;19 int dim,materials_type; 20 20 int numberofelements; 21 21 int numberofvertices; … … 23 23 24 24 /*DataSets: */ 25 Elements *elements = NULL;26 Vertices * vertices= NULL;27 Materials *materials = NULL;25 Elements *elements = NULL; 26 Vertices *vertices = NULL; 27 Materials *materials = NULL; 28 28 29 29 /*Fetch parameters: */ … … 32 32 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); 33 33 iomodel->Constant(&control_analysis,InversionIscontrolEnum); 34 iomodel->Constant(&materials_type,MaterialsEnum); 34 35 35 36 /*Did we already create the elements? : */ … … 44 45 ElementsAndVerticesPartitioning(&iomodel->my_elements,&iomodel->my_vertices,iomodel); 45 46 46 /*Fetch data needed: */ 47 iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum); 47 iomodel->FetchData(2,MeshElementsEnum,MeshElementconnectivityEnum); 48 48 #ifdef _HAVE_3D_ 49 49 if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum); … … 51 51 if(control_analysis)iomodel->FetchData(3,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum); 52 52 53 /*Create elements and materials:*/53 /*Create elements*/ 54 54 for (i=0;i<numberofelements;i++){ 55 55 if(iomodel->my_elements[i]){ … … 60 60 else elements->AddObject(new Penta(i+1,i,i,iomodel,nummodels)); 61 61 #endif 62 63 /*Create and add material property to materials dataset: */64 materials->AddObject(new Matice(i+1,i,iomodel));65 62 } 66 63 } 67 64 65 /*Create materials*/ 66 switch(materials_type){ 67 case MaticeEnum: 68 iomodel->FetchData(2,MaterialsRheologyBEnum,MaterialsRheologyNEnum); 69 for (i=0;i<numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); 70 break; 71 case MatdamageiceEnum: 72 iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum); 73 for (i=0;i<numberofelements;i++) if(iomodel->my_elements[i]) materials->AddObject(new Matdamageice(i+1,i,iomodel)); 74 break; 75 default: 76 _error_("Materials "<<EnumToStringx(materials_type)<<" not supported"); 77 } 78 68 79 /*Free data: */ 69 80 iomodel->DeleteData(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum, -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
r13119 r13129 16 16 void UpdateElementsDiagnosticHoriz(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ 17 17 18 int dim ;18 int dim,materials_type; 19 19 int numberofelements; 20 20 bool ismacayealpattyn; … … 32 32 iomodel->Constant(&control_analysis,InversionIscontrolEnum); 33 33 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); 34 iomodel->Constant(&materials_type,MaterialsEnum); 34 35 35 36 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ … … 61 62 iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); 62 63 iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); 63 iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum);64 64 iomodel->FetchDataToInput(elements,VxEnum); 65 65 iomodel->FetchDataToInput(elements,VyEnum); 66 66 if(materials_type==MatdamageiceEnum){ 67 iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum); 68 } 67 69 if (dim==3){ 68 70 iomodel->FetchDataToInput(elements,MeshElementonbedEnum); -
issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp
r13119 r13129 286 286 else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum; 287 287 else if (strcmp(name,"Matice")==0) return MaticeEnum; 288 else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; 288 289 else if (strcmp(name,"Matpar")==0) return MatparEnum; 289 290 else if (strcmp(name,"Node")==0) return NodeEnum; … … 383 384 else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; 384 385 else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; 385 else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;386 386 else stage=4; 387 387 } 388 388 if(stage==4){ 389 if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; 389 if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; 390 else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; 390 391 else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; 391 392 else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
Note:
See TracChangeset
for help on using the changeset viewer.