Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 13128) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 13129) @@ -282,4 +282,5 @@ MacAyeal3dIceFrontEnum, MaticeEnum, + MatdamageiceEnum, MatparEnum, NodeEnum, Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13128) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13129) @@ -109,4 +109,6 @@ ./classes/objects/Materials/Matice.h\ ./classes/objects/Materials/Matice.cpp\ + ./classes/objects/Materials/Matdamageice.h\ + ./classes/objects/Materials/Matdamageice.cpp\ ./classes/objects/Materials/Matpar.h\ ./classes/objects/Materials/Matpar.cpp\ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13129) @@ -30,5 +30,5 @@ this->nodes=NULL; - this->matice=NULL; + this->material=NULL; this->matpar=NULL; this->verticalneighbors=NULL; @@ -49,5 +49,5 @@ Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels) :PentaRef(nummodels) - ,PentaHook(nummodels,index+1,iomodel) //index+1: matice id, iomodel->numberofelements+1: matpar id + ,PentaHook(nummodels,index+1,iomodel) //index+1: material id, iomodel->numberofelements+1: matpar id { //i is the element index @@ -85,5 +85,5 @@ /*initialize pointers:*/ this->nodes=NULL; - this->matice=NULL; + this->material=NULL; this->matpar=NULL; this->verticalneighbors=NULL; @@ -107,5 +107,5 @@ penta->hnodes=new Hook*[penta->numanalyses]; for(i=0;inumanalyses;i++)penta->hnodes[i]=(Hook*)this->hnodes[i]->copy(); - penta->hmatice=(Hook*)this->hmatice->copy(); + penta->hmaterial=(Hook*)this->hmaterial->copy(); penta->hmatpar=(Hook*)this->hmatpar->copy(); penta->hneighbors=(Hook*)this->hneighbors->copy(); @@ -132,5 +132,5 @@ penta->nodes=xNew(6); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. for(i=0;i<6;i++)penta->nodes[i]=this->nodes[i]; - penta->matice=(Matice*)penta->hmatice->delivers(); + penta->material=(Material*)penta->hmaterial->delivers(); penta->matpar=(Matpar*)penta->hmatpar->delivers(); penta->verticalneighbors=(Penta**)penta->hneighbors->deliverp(); @@ -286,5 +286,5 @@ /*Compute strain rate viscosity and pressure: */ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); pressure_input->GetInputValue(&pressure,gauss); @@ -353,5 +353,5 @@ /*Compute strain rate viscosity and pressure: */ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3d(&viscosity,&epsilon[0]); + material->GetViscosity3d(&viscosity,&epsilon[0]); pressure_input->GetInputValue(&pressure,gauss); @@ -391,5 +391,5 @@ * datasets, using internal ids and offsets hidden in hooks: */ if (this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); - this->hmatice->configure(materialsin); + this->hmaterial->configure(materialsin); this->hmatpar->configure(materialsin); this->hneighbors->configure(elementsin); @@ -398,5 +398,5 @@ if (this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); else this->nodes=NULL; - this->matice=(Matice*)this->hmatice->delivers(); + this->material=(Material*)this->hmaterial->delivers(); this->matpar=(Matpar*)this->hmatpar->delivers(); this->verticalneighbors=(Penta**)this->hneighbors->deliverp(); @@ -419,5 +419,5 @@ /*Checks in debugging {{{*/ - _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs); + _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); /*}}}*/ @@ -490,5 +490,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateKMatrixPrognostic(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Delete Vx and Vy averaged*/ @@ -507,5 +507,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateMassMatrix(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean up and return*/ @@ -522,5 +522,5 @@ /*if debugging mode, check that all pointers exist {{{*/ - _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs); + _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); /*}}}*/ @@ -591,5 +591,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementVector* pe=tria->CreatePVectorPrognostic(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Delete Vx and Vy averaged*/ @@ -609,5 +609,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementVector* pe=tria->CreatePVectorSlope(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean up and return*/ @@ -624,5 +624,5 @@ /*Checks in debugging {{{*/ - _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs); + _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); /*}}}*/ @@ -661,5 +661,5 @@ nodes[4]->DeepEcho(); nodes[5]->DeepEcho(); - matice->DeepEcho(); + material->DeepEcho(); matpar->DeepEcho(); _printLine_(" neighbor ids: " << verticalneighbors[0]->Id() << "-" << verticalneighbors[1]->Id()); @@ -1388,5 +1388,5 @@ original_input=(Input*)penta->inputs->GetInput(enum_type); else if (object_enum==MaterialsEnum) - original_input=(Input*)penta->matice->inputs->GetInput(enum_type); + original_input=(Input*)penta->material->inputs->GetInput(enum_type); else _error_("object " << EnumToStringx(object_enum) << " not supported yet"); @@ -1448,5 +1448,5 @@ this->inputs->AddInput((Input*)depth_averaged_input); else if (object_enum==MaterialsEnum) - this->matice->inputs->AddInput((Input*)depth_averaged_input); + this->material->inputs->AddInput((Input*)depth_averaged_input); else _error_("object " << EnumToStringx(object_enum) << " not supported yet"); @@ -1481,5 +1481,5 @@ num_inputs=1; base_inputs=xNew(num_inputs); - base_inputs[0]=(Input*)matice->inputs->GetInput(enum_type); + base_inputs[0]=(Input*)material->inputs->GetInput(enum_type); } else if (object_type==NodeEnum){ @@ -1515,5 +1515,5 @@ } else if(object_type==MaterialsEnum){ - penta->matice->inputs->AddInput((Input*)copy); + penta->material->inputs->AddInput((Input*)copy); } else if(object_type==NodeEnum){ @@ -1554,6 +1554,6 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyBEnum); - else if (enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyZEnum); + if (enum_type==MaterialsRheologyBbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyBEnum); + else if (enum_type==MaterialsRheologyZbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyZEnum); else input=this->inputs->GetInput(enum_type); //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,5 +1669,5 @@ case MaterialsRheologyBbarEnum: case MaterialsRheologyZbarEnum: - /*Matice will take care of it*/ break; + /*Material will take care of it*/ break; default: _error_("Control " << EnumToStringx(reCast(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); @@ -1965,5 +1965,5 @@ /*update input*/ if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - matice->inputs->AddInput(new PentaP1Input(name,values)); + material->inputs->AddInput(new PentaP1Input(name,values)); } else{ @@ -2637,11 +2637,11 @@ this->SpawnTriaHook(dynamic_cast(tria),&indices[0]); - /*Spawn matice*/ - tria->matice=NULL; - tria->matice=(Matice*)this->matice->copy(); - delete tria->matice->inputs; - tria->matice->inputs=(Inputs*)this->matice->inputs->SpawnTriaInputs(indices); - - /*recover nodes, matice and matpar: */ + /*Spawn material*/ + tria->material=NULL; + tria->material=(Material*)this->material->copy(); + delete tria->material->inputs; + tria->material->inputs=(Inputs*)this->material->inputs->SpawnTriaInputs(indices); + + /*recover nodes, material and matpar: */ tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp(); tria->matpar=(Matpar*)tria->hmatpar->delivers(); @@ -2730,5 +2730,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). S=tria->SurfaceArea(); - delete tria->matice; delete tria; + delete tria->material; delete tria; return S; } @@ -2737,5 +2737,5 @@ tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). S=tria->SurfaceArea(); - delete tria->matice; delete tria; + delete tria->material; delete tria; return S; } @@ -3018,5 +3018,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); GetPhi(&phi, &epsilon[0], viscosity); @@ -3128,5 +3128,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. mass_flux=tria->MassFlux(segment,process_units); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Delete Vx and Vy averaged*/ @@ -3235,8 +3235,8 @@ switch(response_enum){ case MaterialsRheologyBbarEnum: - *presponse=this->matice->GetBbar(); + *presponse=this->material->GetBbar(); break; case MaterialsRheologyZbarEnum: - *presponse=this->matice->GetZbar(); + *presponse=this->material->GetZbar(); break; case VelEnum: @@ -3531,5 +3531,5 @@ ElementMatrix* Ke=tria->CreateKMatrixMelting(); - delete tria->matice; delete tria; + delete tria->material; delete tria; return Ke; } @@ -3830,5 +3830,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); GetPhi(&phi, &epsilon[0], viscosity); @@ -4086,5 +4086,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); GetPhi(&phi, &epsilon[0], viscosity); @@ -4356,5 +4356,5 @@ B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); for(i=0;imatice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); break; case ArrheniusEnum: @@ -4362,7 +4362,7 @@ B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, 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), - matice->GetN()); + material->GetN()); for(i=0;imatice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); break; default: @@ -4434,5 +4434,5 @@ B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0); for(i=0;imatice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); break; case ArrheniusEnum: @@ -4440,7 +4440,7 @@ B_average=Arrhenius((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0, 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), - matice->GetN()); + material->GetN()); for(i=0;imatice->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); + this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B)); break; default: @@ -4468,9 +4468,9 @@ if(enum_type==MaterialsRheologyBbarEnum){ if(!IsOnBed()) return; - input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); + input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); } else if(enum_type==MaterialsRheologyZbarEnum){ if(!IsOnBed()) return; - input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); + input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); } @@ -4491,8 +4491,8 @@ if(enum_type==MaterialsRheologyBbarEnum){ - input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); + input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); } else if(enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); + input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); } else{ @@ -4513,8 +4513,8 @@ if(enum_type==MaterialsRheologyBbarEnum){ - input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); + input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); } else if(enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); + input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); } else{ @@ -4558,20 +4558,25 @@ if (!IsOnBed()) return NULL; - /*Depth Averaging B*/ - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - - /*Depth Averaging Z*/ - this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + /*Depth average some fields*/ + switch(this->material->ObjectEnum()){ + case MaticeEnum: + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + break; + case MatdamageiceEnum: + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + break; + default: + _error_("material "<material->ObjectEnum())<<" not supported"); + } /*Call Tria function*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateKMatrixAdjointMacAyeal(); - delete tria->matice; delete tria; - - /*Delete B averaged*/ - this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); - - /*Delete Z averaged*/ - this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); + delete tria->material; delete tria; + + /*Delete averaged fields*/ + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); /*clean up and return*/ @@ -4619,5 +4624,5 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); + material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; @@ -4689,5 +4694,5 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); + material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; @@ -4754,5 +4759,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementVector* pe=tria->CreatePVectorAdjointHoriz(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean up and return*/ @@ -4768,5 +4773,5 @@ Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). ElementVector* pe=tria->CreatePVectorAdjointHoriz(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean up and return*/ @@ -4782,5 +4787,5 @@ Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). ElementVector* pe=tria->CreatePVectorAdjointStokes(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean up and return*/ @@ -4879,5 +4884,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. tria->GradjDragGradient(gradient,resp,control_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; break; case RheologyBbarAbsGradientEnum: @@ -4885,5 +4890,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. tria->GradjBGradient(gradient,resp,control_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; break; default: @@ -4902,5 +4907,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. tria->GradjDragMacAyeal(gradient,control_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; } /*}}}*/ @@ -5080,8 +5085,8 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). tria->GradjBMacAyeal(gradient,control_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*delete Average B*/ - this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); } /*}}}*/ @@ -5098,8 +5103,8 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); tria->GradjBMacAyeal(gradient,control_index); //We use MacAyeal as an estimate for now - delete tria->matice; delete tria; + delete tria->material; delete tria; /*delete Average B*/ - this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); } /*}}}*/ /*FUNCTION Penta::GradjBbarStokes {{{*/ @@ -5115,8 +5120,8 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); tria->GradjBMacAyeal(gradient,control_index); //We use MacAyeal as an estimate for now - delete tria->matice; delete tria; + delete tria->material; delete tria; /*delete Average B*/ - this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); } /*}}}*/ /*FUNCTION Penta::InputControlUpdate{{{*/ @@ -5136,9 +5141,9 @@ if(control_type[i]==MaterialsRheologyBbarEnum){ if (!IsOnBed()) goto cleanup_and_return; - input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input); + input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input); } else if(control_type[i]==MaterialsRheologyZbarEnum){ if (!IsOnBed()) goto cleanup_and_return; - input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input); + input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input); } else{ @@ -5268,5 +5273,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). J=tria->SurfaceAverageVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5275,5 +5280,5 @@ tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). J=tria->SurfaceAverageVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5305,5 +5310,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). J=tria->SurfaceAbsVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5312,5 +5317,5 @@ tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). J=tria->SurfaceAbsVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5342,5 +5347,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). J=tria->SurfaceLogVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5349,5 +5354,5 @@ tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). J=tria->SurfaceLogVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5381,5 +5386,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). J=tria->SurfaceLogVxVyMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5388,5 +5393,5 @@ tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). J=tria->SurfaceLogVxVyMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5418,5 +5423,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). J=tria->SurfaceRelVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5425,5 +5430,5 @@ tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). J=tria->SurfaceRelVelMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5452,5 +5457,5 @@ tria=(Tria*)SpawnTria(0,1,2); J=tria->ThicknessAbsMisfit(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5467,5 +5472,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria J=tria->DragCoefficientAbsGradient(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5482,5 +5487,5 @@ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria J=tria->RheologyBbarAbsGradient(process_units,weight_index); - delete tria->matice; delete tria; + delete tria->material; delete tria; return J; } @@ -5531,5 +5536,5 @@ if(control_enum==MaterialsRheologyBbarEnum){ - input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); + input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); } else{ @@ -5826,6 +5831,6 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + material->GetViscosity3d(&viscosity, &epsilon[0]); + material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); @@ -5847,5 +5852,5 @@ /*Clean-up and return*/ - delete tria->matice; delete tria; + delete tria->material; delete tria; delete gauss; delete gauss_tria; @@ -6036,5 +6041,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity, &epsilon[0]); + material->GetViscosity3dStokes(&viscosity, &epsilon[0]); D_scalar=2*viscosity*gauss->weight*Jdet; @@ -6063,5 +6068,5 @@ /*Clean-up and return*/ - delete tria->matice; delete tria; + delete tria->material; delete tria; delete gauss; delete gauss_tria; @@ -6143,5 +6148,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); BedNormal(&bed_normal[0],xyz_list_tria); @@ -6335,20 +6340,25 @@ if (!IsOnBed()) return NULL; - /*Depth Averaging B*/ - this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - - /*Depth Averaging Z*/ - this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + /*Depth average some fields*/ + switch(this->material->ObjectEnum()){ + case MaticeEnum: + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + break; + case MatdamageiceEnum: + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + break; + default: + _error_("material "<material->ObjectEnum())<<" not supported"); + } /*Call Tria function*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyeal(); - delete tria->matice; delete tria; - - /*Delete B averaged*/ - this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); - - /*Delete Z averaged*/ - this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); + delete tria->material; delete tria; + + /*Delete averaged fields*/ + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); /*clean up and return*/ @@ -6425,6 +6435,6 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + material->GetViscosity3d(&viscosity, &epsilon[0]); + material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); @@ -6432,5 +6442,5 @@ else if (approximation==MacAyealStokesApproximationEnum){ this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&newviscosity,&epsilons[0]); + material->GetViscosity3dStokes(&newviscosity,&epsilons[0]); } else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); @@ -6451,5 +6461,5 @@ /*Clean up and return*/ - delete tria->matice; + delete tria->material; delete tria; delete gauss_tria; @@ -6469,5 +6479,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean-up and return*/ @@ -6582,5 +6592,5 @@ /*Clean up and return*/ - delete tria->matice; + delete tria->material; delete tria; delete gauss_tria; @@ -6600,5 +6610,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*clean-up and return*/ @@ -6665,6 +6675,6 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + material->GetViscosity3d(&viscosity, &epsilon[0]); + material->GetViscosity3d(&oldviscosity, &oldepsilon[0]); newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); @@ -6829,5 +6839,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); D_scalar=gauss->weight*Jdet; @@ -6900,5 +6910,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); @@ -7088,5 +7098,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); for(i=0;iGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); @@ -7239,5 +7249,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); for(i=0;iGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); @@ -7435,6 +7445,6 @@ rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); - n=matice->GetN(); - B=matice->GetB(); + n=material->GetN(); + B=material->GetB(); Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); @@ -7504,5 +7514,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Clean up and return*/ @@ -7518,5 +7528,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Clean up and return*/ @@ -7636,5 +7646,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); for(i=0;iInputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); - - /*Depth Averaging Z*/ - this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + /*Depth average some fields*/ + switch(this->material->ObjectEnum()){ + case MaticeEnum: + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + break; + case MatdamageiceEnum: + this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum); + this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum); + break; + default: + _error_("material "<material->ObjectEnum())<<" not supported"); + } /*Call Tria function*/ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal(); - delete tria->matice; delete tria; - - /*Delete B averaged*/ - this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum); - - /*Delete Z averaged*/ - this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum); + delete tria->material; delete tria; + + /*Delete averaged inputs*/ + this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum); + this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum); /*clean up and return*/ @@ -7955,5 +7970,5 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); + material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; @@ -8022,5 +8037,5 @@ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); + material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3]; eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4]; @@ -8264,6 +8279,6 @@ /*Get A*/ - _assert_(matice->GetN()==3.0); - A=matice->GetA(); + _assert_(material->GetN()==3.0); + A=material->GetA(); /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/ @@ -9085,5 +9100,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementMatrix* Ke=tria->CreateKMatrixBalancethickness(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Delete Vx and Vy averaged*/ @@ -9107,5 +9122,5 @@ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. ElementVector* pe=tria->CreatePVectorBalancethickness(); - delete tria->matice; delete tria; + delete tria->material; delete tria; /*Delete Vx and Vy averaged*/ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 13129) @@ -16,5 +16,5 @@ class IoModel; class Node; -class Matice; +class Material; class Matpar; class Tria; @@ -35,5 +35,5 @@ Node **nodes; // 6 nodes - Matice *matice; // 1 material ice + Material *material; // 1 material ice Matpar *matpar; // 1 material parameter Penta **verticalneighbors; // 2 neighbors: first one under, second one above Index: /issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.cpp (revision 13129) @@ -25,5 +25,5 @@ numanalyses=UNDEF; this->hnodes=NULL; - this->hmatice=NULL; + this->hmaterial=NULL; this->hmatpar=NULL; this->hneighbors=NULL; @@ -39,11 +39,11 @@ } delete [] this->hnodes; - delete hmatice; + delete hmaterial; delete hmatpar; delete hneighbors; } /*}}}*/ -/*FUNCTION PentaHook::PentaHook(int in_numanalyses,int matice_id, int matpar_id){{{*/ -PentaHook::PentaHook(int in_numanalyses,int matice_id, IoModel* iomodel){ +/*FUNCTION PentaHook::PentaHook(int in_numanalyses,int material_id, int matpar_id){{{*/ +PentaHook::PentaHook(int in_numanalyses,int material_id, IoModel* iomodel){ /*intermediary: */ @@ -55,5 +55,5 @@ this->numanalyses=in_numanalyses; this->hnodes=new Hook*[in_numanalyses]; - this->hmatice=new Hook(&matice_id,1); + this->hmaterial=new Hook(&material_id,1); this->hmatpar=new Hook(&matpar_id,1); this->hneighbors=NULL; @@ -98,6 +98,6 @@ } } - // do not spawn hmatice. matice will be taken care of by Penta - triahook->hmatice=NULL; + // do not spawn hmaterial. material will be taken care of by Penta + triahook->hmaterial=NULL; triahook->hmatpar=(Hook*)this->hmatpar->copy(); } Index: /issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/PentaHook.h (revision 13129) @@ -15,5 +15,5 @@ int numanalyses; //number of analysis types Hook** hnodes; // 6 nodes for each analysis type - Hook* hmatice; // 1 ice material + Hook* hmaterial; // 1 ice material Hook* hmatpar; // 1 material parameter Hook* hneighbors; // 2 elements, first down, second up @@ -21,5 +21,5 @@ /*FUNCTION constructors, destructors {{{*/ PentaHook(); - PentaHook(int in_numanalyses,int matice_id, IoModel* iomodel); + PentaHook(int in_numanalyses,int material_id, IoModel* iomodel); ~PentaHook(); void SetHookNodes(int* node_ids,int analysis_counter); Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13129) @@ -29,5 +29,5 @@ this->nodes=NULL; - this->matice=NULL; + this->material=NULL; this->matpar=NULL; for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF; @@ -61,5 +61,5 @@ /*initialize pointers:*/ this->nodes=NULL; - this->matice=NULL; + this->material=NULL; this->matpar=NULL; @@ -89,5 +89,5 @@ tria->hnodes=new Hook*[tria->numanalyses]; for(i=0;inumanalyses;i++)tria->hnodes[i]=(Hook*)this->hnodes[i]->copy(); - tria->hmatice=(Hook*)this->hmatice->copy(); + tria->hmaterial=(Hook*)this->hmaterial->copy(); tria->hmatpar=(Hook*)this->hmatpar->copy(); @@ -113,5 +113,5 @@ tria->nodes=xNew(3); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. for(i=0;i<3;i++)tria->nodes[i]=this->nodes[i]; - tria->matice=(Matice*)tria->hmatice->delivers(); + tria->material=(Material*)tria->hmaterial->delivers(); tria->matpar=(Matpar*)tria->hmatpar->delivers(); @@ -173,5 +173,5 @@ /*Checks in debugging mode{{{*/ - _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); /*}}}*/ @@ -526,5 +526,5 @@ /*asserts: {{{*/ /*if debugging mode, check that all pointers exist*/ - _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); /*}}}*/ @@ -748,5 +748,5 @@ /*Checks in debugging {{{*/ - _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs); /*}}}*/ @@ -812,5 +812,5 @@ /*Compute strain rate viscosity and pressure: */ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosity2d(&viscosity,&epsilon[0]); + material->GetViscosity2d(&viscosity,&epsilon[0]); pressure_input->GetInputValue(&pressure,gauss); @@ -846,5 +846,5 @@ * datasets, using internal ids and offsets hidden in hooks: */ if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); - this->hmatice->configure(materialsin); + this->hmaterial->configure(materialsin); this->hmatpar->configure(materialsin); @@ -852,5 +852,5 @@ if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); else this->nodes=NULL; - this->matice=(Matice*)this->hmatice->delivers(); + this->material=(Material*)this->hmaterial->delivers(); this->matpar=(Matpar*)this->hmatpar->delivers(); @@ -875,6 +875,6 @@ else _printLine_("nodes = NULL"); - if (matice) matice->DeepEcho(); - else _printLine_("matice = NULL"); + if (material) material->DeepEcho(); + else _printLine_("material = NULL"); if (matpar) matpar->DeepEcho(); @@ -987,6 +987,6 @@ else _printLine_("nodes = NULL"); - if (matice) matice->Echo(); - else _printLine_("matice = NULL"); + if (material) material->Echo(); + else _printLine_("material = NULL"); if (matpar) matpar->Echo(); @@ -1350,5 +1350,5 @@ oldinput=(Input*)this->inputs->GetInput(enum_type); else if (object_enum==MaterialsEnum) - oldinput=(Input*)this->matice->inputs->GetInput(enum_type); + oldinput=(Input*)this->material->inputs->GetInput(enum_type); else _error_("object " << EnumToStringx(object_enum) << " not supported yet"); @@ -1363,5 +1363,5 @@ this->inputs->AddInput((Input*)newinput); else if (object_enum==MaterialsEnum) - this->matice->inputs->AddInput((Input*)newinput); + this->material->inputs->AddInput((Input*)newinput); else _error_("object " << EnumToStringx(object_enum) << " not supported yet"); @@ -1396,5 +1396,5 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type); + if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->material->inputs->GetInput(enum_type); else input=this->inputs->GetInput(enum_type); //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs"); @@ -1503,5 +1503,5 @@ case MaterialsRheologyBbarEnum: case MaterialsRheologyZbarEnum: - /*Matice will take care of it*/ break; + /*Material will take care of it*/ break; default: _error_("Control " << EnumToStringx(reCast(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet"); @@ -1692,5 +1692,5 @@ /*update input*/ if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ - matice->inputs->AddInput(new TriaP1Input(name,values)); + material->inputs->AddInput(new TriaP1Input(name,values)); } else{ @@ -2792,8 +2792,8 @@ switch(response_enum){ case MaterialsRheologyBbarEnum: - *presponse=this->matice->GetBbar(); + *presponse=this->material->GetBbar(); break; case MaterialsRheologyZbarEnum: - *presponse=this->matice->GetZbar(); + *presponse=this->material->GetZbar(); break; case VelEnum:{ @@ -2909,6 +2909,6 @@ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); - matice->GetViscosity2d(&viscosity, &epsilon[0]); - matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); + material->GetViscosity2d(&viscosity, &epsilon[0]); + material->GetViscosity2d(&oldviscosity, &oldepsilon[0]); thickness_input->GetInputValue(&thickness, gauss); @@ -3092,6 +3092,6 @@ rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); - n=matice->GetN(); - B=matice->GetBbar(); + n=material->GetN(); + B=material->GetBbar(); Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input); Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input); @@ -3164,5 +3164,5 @@ thickness_input->GetInputValue(&thickness, gauss); this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); + material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; @@ -3409,5 +3409,5 @@ if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input); + input=(Input*)material->inputs->GetInput(control_type[i]); _assert_(input); } else{ @@ -3436,5 +3436,5 @@ if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(enum_type); + input=(Input*)material->inputs->GetInput(enum_type); } else{ @@ -3454,5 +3454,5 @@ if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(enum_type); + input=(Input*)material->inputs->GetInput(enum_type); } else{ @@ -3473,5 +3473,5 @@ if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(enum_type); + input=(Input*)material->inputs->GetInput(enum_type); } else{ @@ -3567,5 +3567,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GradientIndexing(&doflist1[0],control_index); - Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); @@ -3607,5 +3607,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GradientIndexing(&doflist1[0],control_index); - Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + Input* rheologyz_input=material->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); @@ -3656,5 +3656,5 @@ Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); /* Start looping on the number of gaussian points: */ @@ -3672,5 +3672,5 @@ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); + material->GetViscosityComplement(&viscosity_complement,&epsilon[0]); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); @@ -3713,5 +3713,5 @@ Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + Input* rheologyz_input=material->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); /* Start looping on the number of gaussian points: */ @@ -3729,5 +3729,5 @@ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]); + material->GetViscosityZComplement(&viscosity_complement,&epsilon[0]); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); @@ -4013,5 +4013,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); - Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); + Input* rheologyb_input=material->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input); /* Start looping on the number of gaussian points: */ @@ -5110,5 +5110,5 @@ thickness_input->GetInputValue(&thickness, gauss); this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); + material->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]); eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2]; eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1]; @@ -5216,5 +5216,5 @@ /*Get input (either in element or material)*/ if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); + input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); } else{ @@ -5251,5 +5251,5 @@ if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){ - input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input); + input=(Input*)material->inputs->GetInput(control_enum); _assert_(input); } else{ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h (revision 13129) @@ -15,5 +15,5 @@ class IoModel; class Node; -class Matice; +class Material; class Matpar; class ElementMatrix; @@ -33,5 +33,5 @@ Node **nodes; // 3 nodes - Matice *matice; // 1 material ice + Material *material; // 1 material ice Matpar *matpar; // 1 material parameter int horizontalneighborsids[3]; Index: /issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.cpp (revision 13129) @@ -25,5 +25,5 @@ numanalyses=UNDEF; this->hnodes=NULL; - this->hmatice=NULL; + this->hmaterial=NULL; this->hmatpar=NULL; } @@ -37,11 +37,11 @@ } delete [] this->hnodes; - delete hmatice; + delete hmaterial; delete hmatpar; } /*}}}*/ -/*FUNCTION TriaHook::TriaHook(int in_numanalyses,int matice_id, int matpar_id){{{*/ -TriaHook::TriaHook(int in_numanalyses,int matice_id, IoModel* iomodel){ +/*FUNCTION TriaHook::TriaHook(int in_numanalyses,int material_id, int matpar_id){{{*/ +TriaHook::TriaHook(int in_numanalyses,int material_id, IoModel* iomodel){ /*intermediary: */ @@ -53,5 +53,5 @@ this->numanalyses=in_numanalyses; this->hnodes= new Hook*[in_numanalyses]; - this->hmatice=new Hook(&matice_id,1); + this->hmaterial=new Hook(&material_id,1); this->hmatpar=new Hook(&matpar_id,1); Index: /issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/TriaHook.h (revision 13129) @@ -14,5 +14,5 @@ int numanalyses; //number of analysis types Hook** hnodes; // 3 nodes for each analysis type - Hook* hmatice; // 1 ice material + Hook* hmaterial; // 1 ice material Hook* hmatpar; // 1 material parameter @@ -20,5 +20,5 @@ /*FUNCTION constructors, destructors {{{*/ TriaHook(); - TriaHook(int in_numanalyses,int matice_id, IoModel* iomodel); + TriaHook(int in_numanalyses,int material_id, IoModel* iomodel); ~TriaHook(); void SetHookNodes(int* node_ids,int analysis_counter); Index: /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Loads/Icefront.cpp (revision 13129) @@ -498,5 +498,5 @@ /*clean-up and return*/ - delete tria->matice; + delete tria->material; delete tria; delete icefront; Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Material.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Material.h (revision 13129) @@ -17,10 +17,25 @@ public: + Inputs* inputs; virtual ~Material(){}; + /*WARNING: input should not be public but it is an easy way to update B from T (using UpdateFromSolution) from Pentas*/ /*Numerics*/ - virtual void InputDuplicate(int original_enum,int new_enum)=0; - virtual void Configure(Elements* elements)=0; - virtual void GetVectorFromInputs(Vector* vector,int input_enum)=0; + virtual void InputDuplicate(int original_enum,int new_enum)=0; + virtual void Configure(Elements* elements)=0; + virtual void GetVectorFromInputs(Vector* vector,int input_enum)=0; + virtual void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0; + virtual void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0; + virtual void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0; + virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; + virtual void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0; + virtual void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; + virtual void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0; + virtual IssmDouble GetA()=0; + virtual IssmDouble GetB()=0; + virtual IssmDouble GetBbar()=0; + virtual IssmDouble GetN()=0; + virtual IssmDouble GetZ()=0; + virtual IssmDouble GetZbar()=0; }; Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp (revision 13129) @@ -172,24 +172,4 @@ } /*}}}*/ -/*FUNCTION Matice::GetZ {{{*/ -IssmDouble Matice::GetZ(){ - - /*Output*/ - IssmDouble Z; - - inputs->GetInputAverage(&Z,MaterialsRheologyZEnum); - return Z; -} -/*}}}*/ -/*FUNCTION Matice::GetZbar {{{*/ -IssmDouble Matice::GetZbar(){ - - /*Output*/ - IssmDouble Zbar; - - inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum); - return Zbar; -} -/*}}}*/ /*FUNCTION Matice::GetVectorFromInputs{{{*/ void Matice::GetVectorFromInputs(Vector* vector,int input_enum){ @@ -227,5 +207,5 @@ void Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ /*From a string tensor and a material object, return viscosity, using Glen's flow law. - Z * B + B viscosity= ------------------------------------------------------------------- 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -246,11 +226,9 @@ /*Intermediary: */ IssmDouble A,e; - IssmDouble Btmp,B,n,Z; + IssmDouble B,n; /*Get B and n*/ - Btmp=GetBbar(); - Z=GetZbar(); + B=GetBbar(); n=GetN(); - B=Z*Btmp; if (n==1){ @@ -314,10 +292,9 @@ /*Intermediaries: */ IssmDouble A,e; - IssmDouble B,n,Z; - - /*Get B, Z and n*/ + IssmDouble B,n; + + /*Get B and n*/ + B=GetB(); n=GetN(); - Z=GetZ(); - B=Z*GetB(); if (n==1){ @@ -385,12 +362,11 @@ /*Intermediaries: */ IssmDouble A,e; - IssmDouble B,n,Z; + IssmDouble B,n; IssmDouble eps0; /*Get B and n*/ eps0=pow((IssmDouble)10,(IssmDouble)-27); + B=GetB(); n=GetN(); - Z=GetZ(); - B=Z*GetB(); if (n==1){ @@ -476,60 +452,4 @@ viscosity_complement=1/(2*pow(A,e)); - } - } - else{ - viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); - } - - /*Checks in debugging mode*/ - _assert_(B>0); - _assert_(n>0); - _assert_(viscosity_complement>0); - - /*Return: */ - *pviscosity_complement=viscosity_complement; -} -/*}}}*/ -/*FUNCTION Matice::GetViscosityZComplement {{{*/ -void Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ - /*Return viscosity derivative for control method d(mu)/dZ: - * - * B - * dviscosity= ------------------------------------------------------------------- - * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] - * - * If epsilon is NULL, it means this is the first time Gradjb is being run, and we - * return mu20, initial viscosity. - */ - - /*output: */ - IssmDouble viscosity_complement; - - /*input strain rate: */ - IssmDouble exx,eyy,exy; - - /*Intermediary value A and exponent e: */ - IssmDouble A,e; - IssmDouble B,n; - - /*Get B and n*/ - B=GetBbar(); - n=GetN(); - - if(epsilon){ - exx=*(epsilon+0); - eyy=*(epsilon+1); - exy=*(epsilon+2); - - /*Build viscosity: mu2=B/(2*A^e) */ - A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; - if(A==0){ - /*Maximum viscosity_complement for 0 shear areas: */ - viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); - } - else{ - e=(n-1)/(2*n); - - viscosity_complement=B/(2*pow(A,e)); } } @@ -781,10 +701,4 @@ } - /*Get Z*/ - if (iomodel->Data(MaterialsRheologyZEnum)) { - for(i=0;iData(MaterialsRheologyZEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); - } - /*Control Inputs*/ #ifdef _HAVE_CONTROL_ @@ -801,14 +715,4 @@ } break; - case MaterialsRheologyZbarEnum: - if (iomodel->Data(MaterialsRheologyZEnum)){ - _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); - for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; - for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); - } - break; - } } @@ -837,10 +741,4 @@ for(i=0;iData(MaterialsRheologyNEnum)[index]; this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); - } - - /*Get Z*/ - if (iomodel->Data(MaterialsRheologyZEnum)) { - for(i=0;iData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; - this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs)); } @@ -859,13 +757,4 @@ } break; - case MaterialsRheologyZbarEnum: - if (iomodel->Data(MaterialsRheologyZEnum)){ - _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); - for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; - for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; - this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); - } - break; } } @@ -886,7 +775,5 @@ name==MaterialsRheologyBEnum || name==MaterialsRheologyBbarEnum || - name==MaterialsRheologyNEnum || - name==MaterialsRheologyZEnum || - name==MaterialsRheologyZbarEnum + name==MaterialsRheologyNEnum ){ return true; Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h (revision 13129) @@ -15,14 +15,8 @@ private: - /*Id*/ int mid; - - /*hooks: */ Hook* helement; public: - /*WARNING: input should not be public but it is an easy way to update B from T (using UpdateFromSolution) from Pentas*/ - Inputs* inputs; - /*Matice constructors, destructors: {{{*/ Matice(); @@ -56,21 +50,19 @@ void Configure(Elements* elements); void GetVectorFromInputs(Vector* vector,int input_enum); - /*}}}*/ - /*Matice Numerics: {{{*/ - void SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); - void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); - void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); - void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); - void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); - void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); + void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); + void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); + void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); + void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityZComplement(IssmDouble*, IssmDouble*){_error_("not supported");}; + void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); IssmDouble GetA(); IssmDouble GetB(); IssmDouble GetBbar(); + IssmDouble GetZ(){_error_("not supported");}; + IssmDouble GetZbar(){_error_("not supported");}; IssmDouble GetN(); - IssmDouble GetZ(); - IssmDouble GetZbar(); - bool IsInput(int name); + bool IsInput(int name); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp (revision 13129) @@ -45,4 +45,5 @@ iomodel->Constant(&this->hydro_p,HydrologyPEnum); iomodel->Constant(&this->hydro_q,HydrologyQEnum); + this->inputs=NULL; /*not used here*/ } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h (revision 13129) @@ -68,4 +68,17 @@ void Configure(Elements* elements); void GetVectorFromInputs(Vector* vector,int input_enum){return;} + void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");}; + void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; + void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");}; + IssmDouble GetA(){_error_("not supported");}; + IssmDouble GetB(){_error_("not supported");}; + IssmDouble GetBbar(){_error_("not supported");}; + IssmDouble GetN(){_error_("not supported");}; + IssmDouble GetZ(){_error_("not supported");}; + IssmDouble GetZbar(){_error_("not supported");}; /*}}}*/ /*Numerics: {{{*/ Index: /issm/trunk-jpl/src/c/classes/objects/Params/BoolParam.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Params/BoolParam.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/Params/BoolParam.h (revision 13129) @@ -43,31 +43,31 @@ int InstanceEnum(){return enum_type;} void GetParameterValue(bool* pbool){*pbool=value;} - void GetParameterValue(int* pinteger){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return an integer");} - void GetParameterValue(int** pintarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return an array of integers");} - void GetParameterValue(int** pintarray,int* pM,int* pN){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return an array of integers");} - void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a IssmDouble");} + void GetParameterValue(int* pinteger){_error_("Param "<< EnumToStringx(enum_type) << " cannot return an integer");} + void GetParameterValue(int** pintarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return an array of integers");} + void GetParameterValue(int** pintarray,int* pM,int* pN){_error_("Param "<< EnumToStringx(enum_type) << " cannot return an array of integers");} + void GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");} void GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");} - void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a string");} - void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a string array");} - void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a IssmDouble array");} - void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, int* pN){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a IssmDouble array");} - void GetParameterValue(IssmDouble*** parray, int* pM,int** pmdims, int** pndims){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a matrix array");} - void GetParameterValue(Vector** pvec){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a Vec");} - void GetParameterValue(Matrix** pmat){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a Mat");} - void GetParameterValue(FILE** pfid){_error_("Param "<< EnumToStringx(enum_type) << ") cannot return a FILE");} + void GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");} + void GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble array");} + void GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, int* pN){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble array");} + void GetParameterValue(IssmDouble*** parray, int* pM,int** pmdims, int** pndims){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a matrix array");} + void GetParameterValue(Vector** pvec){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a Vec");} + void GetParameterValue(Matrix** pmat){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a Mat");} + void GetParameterValue(FILE** pfid){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a FILE");} void SetValue(bool boolean){this->value=boolean;} - void SetValue(int integer){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold an int");} - void SetValue(IssmDouble scalar){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold an IssmPDouble");} - void SetValue(char* string){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a string");} - void SetValue(char** stringarray,int M){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a string array");} - void SetValue(IssmDouble* IssmDoublearray,int M){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a IssmDouble array");} - void SetValue(IssmDouble* pIssmDoublearray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a IssmDouble array");} - void SetValue(int* intarray,int M){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a int array");} - void SetValue(int* pintarray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a int array");} - void SetValue(Vector* vec){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a Vec");} - void SetValue(Matrix* mat){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a Mat");} - void SetValue(FILE* fid){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold a FILE");} - void SetValue(IssmDouble** array, int M, int* mdim_array, int* ndim_array){_error_("Param "<< EnumToStringx(enum_type) << ") cannot hold an array of matrices");} + void SetValue(int integer){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an int");} + void SetValue(IssmDouble scalar){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an IssmPDouble");} + void SetValue(char* string){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a string");} + void SetValue(char** stringarray,int M){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a string array");} + void SetValue(IssmDouble* IssmDoublearray,int M){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a IssmDouble array");} + void SetValue(IssmDouble* pIssmDoublearray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a IssmDouble array");} + void SetValue(int* intarray,int M){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a int array");} + void SetValue(int* pintarray,int M,int N){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a int array");} + void SetValue(Vector* vec){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a Vec");} + void SetValue(Matrix* mat){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a Mat");} + void SetValue(FILE* fid){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold a FILE");} + void SetValue(IssmDouble** array, int M, int* mdim_array, int* ndim_array){_error_("Param "<< EnumToStringx(enum_type) << " cannot hold an array of matrices");} void UnitConversion(int direction_enum); Index: /issm/trunk-jpl/src/c/classes/objects/objects.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/objects.h (revision 13128) +++ /issm/trunk-jpl/src/c/classes/objects/objects.h (revision 13129) @@ -108,4 +108,5 @@ #include "./Materials/Material.h" #include "./Materials/Matice.h" +#include "./Materials/Matdamageice.h" #include "./Materials/Matpar.h" Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 13129) @@ -279,4 +279,5 @@ case MacAyeal3dIceFrontEnum : return "MacAyeal3dIceFront"; case MaticeEnum : return "Matice"; + case MatdamageiceEnum : return "Matdamageice"; case MatparEnum : return "Matpar"; case NodeEnum : return "Node"; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 13129) @@ -17,5 +17,5 @@ /*Intermediary*/ int i,j,k,n; - int dim; + int dim,materials_type; int numberofelements; int numberofvertices; @@ -23,7 +23,7 @@ /*DataSets: */ - Elements* elements = NULL; - Vertices* vertices = NULL; - Materials* materials = NULL; + Elements *elements = NULL; + Vertices *vertices = NULL; + Materials *materials = NULL; /*Fetch parameters: */ @@ -32,4 +32,5 @@ iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); iomodel->Constant(&control_analysis,InversionIscontrolEnum); + iomodel->Constant(&materials_type,MaterialsEnum); /*Did we already create the elements? : */ @@ -44,6 +45,5 @@ ElementsAndVerticesPartitioning(&iomodel->my_elements,&iomodel->my_vertices,iomodel); - /*Fetch data needed: */ - iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum); + iomodel->FetchData(2,MeshElementsEnum,MeshElementconnectivityEnum); #ifdef _HAVE_3D_ if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum); @@ -51,5 +51,5 @@ if(control_analysis)iomodel->FetchData(3,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum); - /*Create elements and materials: */ + /*Create elements*/ for (i=0;imy_elements[i]){ @@ -60,10 +60,21 @@ else elements->AddObject(new Penta(i+1,i,i,iomodel,nummodels)); #endif - - /*Create and add material property to materials dataset: */ - materials->AddObject(new Matice(i+1,i,iomodel)); } } + /*Create materials*/ + switch(materials_type){ + case MaticeEnum: + iomodel->FetchData(2,MaterialsRheologyBEnum,MaterialsRheologyNEnum); + for (i=0;imy_elements[i]) materials->AddObject(new Matice(i+1,i,iomodel)); + break; + case MatdamageiceEnum: + iomodel->FetchData(3,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum); + for (i=0;imy_elements[i]) materials->AddObject(new Matdamageice(i+1,i,iomodel)); + break; + default: + _error_("Materials "<DeleteData(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum, Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 13129) @@ -16,5 +16,5 @@ void UpdateElementsDiagnosticHoriz(Elements* elements, IoModel* iomodel,int analysis_counter,int analysis_type){ - int dim; + int dim,materials_type; int numberofelements; bool ismacayealpattyn; @@ -32,4 +32,5 @@ iomodel->Constant(&control_analysis,InversionIscontrolEnum); iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); + iomodel->Constant(&materials_type,MaterialsEnum); /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ @@ -61,8 +62,9 @@ iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); - iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum); iomodel->FetchDataToInput(elements,VxEnum); iomodel->FetchDataToInput(elements,VyEnum); - + if(materials_type==MatdamageiceEnum){ + iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum); + } if (dim==3){ iomodel->FetchDataToInput(elements,MeshElementonbedEnum); Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 13128) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 13129) @@ -286,4 +286,5 @@ else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum; else if (strcmp(name,"Matice")==0) return MaticeEnum; + else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matpar")==0) return MatparEnum; else if (strcmp(name,"Node")==0) return NodeEnum; @@ -383,9 +384,9 @@ else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; - else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; + if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; + else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;