Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15400) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15401) @@ -2178,4 +2178,5 @@ name==FrictionCoefficientEnum || name==GLlevelsetEnum || + name==IcelevelsetEnum || name==GradientEnum || name==OldGradientEnum || Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15400) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15401) @@ -1733,4 +1733,5 @@ name==BedEnum || name==GLlevelsetEnum || + name==IcelevelsetEnum || name==SurfaceSlopeXEnum || name==SurfaceSlopeYEnum || @@ -2942,4 +2943,18 @@ ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){ + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorDiagnosticMacAyealDrivingStress(); + ElementVector* pe2=CreatePVectorDiagnosticMacAyealFront(); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorDiagnosticMacAyealDrivingStress {{{*/ +ElementVector* Tria::CreatePVectorDiagnosticMacAyealDrivingStress(){ + /*Intermediaries */ int i,j; @@ -2948,4 +2963,5 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble slope[2]; + IssmDouble icefrontlevel[3]; /*Fetch number of nodes and dof for this finite element*/ @@ -2963,4 +2979,5 @@ Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input); + GetInputListOnVertices(&icefrontlevel[0],BedEnum); /* Start looping on the number of gaussian points: */ @@ -2981,4 +2998,41 @@ pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i]; } + } + } + + /*Transform coordinate system*/ + TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum); + + /*Clean up and return*/ + delete gauss; + xDelete(basis); + return pe; +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorDiagnosticMacAyealFront {{{*/ +ElementVector* Tria::CreatePVectorDiagnosticMacAyealFront(){ + + /*Intermediaries */ + int i,j; + IssmDouble ls[3]; + IssmDouble xyz_list[NUMVERTICES][3]; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + int numdof = numnodes*NDOF2; + Icefront *icefront=NULL; + + return NULL; + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + GetInputListOnVertices(&ls[0],IcelevelsetEnum); + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + GaussTria* gauss = new GaussTria(2); + IssmDouble* basis = xNew(numnodes); + + /*Create Ice Front if necessary*/ + if(ls[0]>0. || ls[1]>0. || ls[2]>0.){ + if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){ + //icefront=new Icefront("2d",inputs,matpar,MacAyealApproximationEnum,analysis_type); } } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15400) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15401) @@ -224,4 +224,6 @@ ElementMatrix* CreateKMatrixDiagnosticHutter(void); ElementVector* CreatePVectorDiagnosticMacAyeal(void); + ElementVector* CreatePVectorDiagnosticMacAyealDrivingStress(void); + ElementVector* CreatePVectorDiagnosticMacAyealFront(void); ElementVector* CreatePVectorDiagnosticHutter(void); ElementMatrix* CreateJacobianDiagnosticMacayeal(void); Index: /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp (revision 15400) +++ /issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp (revision 15401) @@ -122,4 +122,87 @@ this->element = NULL; this->matpar = NULL; +} + +/*}}}*/ +/*FUNCTION Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type) {{{*/ +Icefront::Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in,int in_icefront_type, int in_analysis_type){ + + int segment_width; + int element; + int numnodes; + int numvertices; + int dim; + int numberofelements; + + /*icefront constructor data: */ + int icefront_eid; + int icefront_mparid; + int icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size + int icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size + +// /*find parameters: */ +// iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); +// + /*First, retrieve element index and element type: */ + if(strcmp(element_type_in,"2d")==0){ + segment_width=4; + } + else{ + segment_width=6; + } +// element=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill) +// +// /*Build ids for hook constructors: */ +// icefront_eid=reCast( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing +// icefront_mparid=numberofelements+1; //matlab indexing +// + if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){ +// icefront_node_ids[0]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); +// icefront_node_ids[1]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); +// icefront_vertex_ids[0]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); +// icefront_vertex_ids[1]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); + } + else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){ +// icefront_node_ids[0]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); +// icefront_node_ids[1]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); +// icefront_node_ids[2]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2)); +// icefront_node_ids[3]=iomodel->nodecounter+reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3)); +// icefront_vertex_ids[0]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0)); +// icefront_vertex_ids[1]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1)); +// icefront_vertex_ids[2]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2)); +// icefront_vertex_ids[3]=reCast(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3)); + } + else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!"); + + if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){ + numnodes=4; + numvertices=4; + } + else{ + numnodes=2; + numvertices=2; + } + + /*Ok, we have everything to build the object: */ + this->id=1; + this->analysis_type=in_analysis_type; + + /*Hooks: */ + this->hnodes=new Hook(icefront_node_ids,numnodes); + this->hvertices=new Hook(icefront_vertex_ids,numvertices); + this->helement=new Hook(&icefront_eid,1); + this->hmatpar=new Hook(&icefront_mparid,1); + + //intialize and add as many inputs per element as requested: + this->inputs=inputs_in; + this->inputs->AddInput(new IntInput(FillEnum,1)); //We always consider we have water, if above sea level, only air will be applied + this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type)); + + //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this. + this->parameters = NULL; + this->nodes = NULL; + this->vertices = NULL; + this->element = NULL; + this->matpar = matpar_in; } Index: /issm/trunk-jpl/src/c/classes/Loads/Icefront.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Icefront.h (revision 15400) +++ /issm/trunk-jpl/src/c/classes/Loads/Icefront.h (revision 15401) @@ -45,4 +45,5 @@ Icefront(); Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int analysis_type); + Icefront(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int icefront_type, int in_analysis_type); ~Icefront(); /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 15400) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 15401) @@ -92,5 +92,5 @@ /*Fetch data:*/ - iomodel->FetchData(6,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum); + iomodel->FetchData(7,MeshElementsEnum,MeshXEnum,MeshYEnum,MeshZEnum,BedEnum,ThicknessEnum,IcelevelsetEnum); CreateNumberNodeToElementConnectivity(iomodel); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15400) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 15401) @@ -144,4 +144,5 @@ MaskVertexongroundediceEnum, MaskVertexonwaterEnum, + IcelevelsetEnum, MaterialsBetaEnum, MaterialsHeatcapacityEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15400) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 15401) @@ -152,4 +152,5 @@ case MaskVertexongroundediceEnum : return "MaskVertexongroundedice"; case MaskVertexonwaterEnum : return "MaskVertexonwater"; + case IcelevelsetEnum : return "Icelevelset"; case MaterialsBetaEnum : return "MaterialsBeta"; case MaterialsHeatcapacityEnum : return "MaterialsHeatcapacity"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15400) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 15401) @@ -155,4 +155,5 @@ else if (strcmp(name,"MaskVertexongroundedice")==0) return MaskVertexongroundediceEnum; else if (strcmp(name,"MaskVertexonwater")==0) return MaskVertexonwaterEnum; + else if (strcmp(name,"Icelevelset")==0) return IcelevelsetEnum; else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum; else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; - else if (strcmp(name,"Thickness")==0) return ThicknessEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum; + if (strcmp(name,"Thickness")==0) return ThicknessEnum; + else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum; else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum; else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum; else if (strcmp(name,"RiftfrontType")==0) return RiftfrontTypeEnum; - else if (strcmp(name,"Segment")==0) return SegmentEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; + if (strcmp(name,"Segment")==0) return SegmentEnum; + else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum; else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum; - else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; + if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; + else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else if (strcmp(name,"J")==0) return JEnum; Index: /issm/trunk-jpl/src/m/classes/mask.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.m (revision 15400) +++ /issm/trunk-jpl/src/m/classes/mask.m (revision 15401) @@ -54,4 +54,8 @@ WriteData(fid,'object',obj,'fieldname','vertexongroundedice','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','vertexonwater','format','DoubleMat','mattype',1); + icelevelset=ones(md.mesh.numberofvertices,1); + pos=md.diagnostic.icefront(:,1:end-1); + icelevelset(pos(:))=0; + WriteData(fid,'data',icelevelset,'format','DoubleMat','mattype',1,'enum',IcelevelsetEnum()); end % }}} end Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 15400) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 15401) @@ -1899,4 +1899,18 @@ return StringToEnum('MaskVertexonwater')[0] +def IcelevelsetEnum(): + """ + ICELEVELSETENUM - Enum of Icelevelset + + WARNING: DO NOT MODIFY THIS FILE + this file has been automatically generated by src/c/shared/Enum/Synchronize.sh + Please read src/c/shared/Enum/README for more information + + Usage: + macro=IcelevelsetEnum() + """ + + return StringToEnum('Icelevelset')[0] + def MaterialsBetaEnum(): """ @@ -7861,4 +7875,4 @@ """ - return 560 - + return 561 + Index: /issm/trunk-jpl/src/m/enum/IcelevelsetEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/IcelevelsetEnum.m (revision 15401) +++ /issm/trunk-jpl/src/m/enum/IcelevelsetEnum.m (revision 15401) @@ -0,0 +1,11 @@ +function macro=IcelevelsetEnum() +%ICELEVELSETENUM - Enum of Icelevelset +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=IcelevelsetEnum() + +macro=StringToEnum('Icelevelset'); Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 15400) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 15401) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=560; +macro=561;