Index: /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 3946) +++ /issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp (revision 3947) @@ -161,4 +161,5 @@ case NodeOnIceShelfEnum : return "NodeOnIceShelf"; case NodeOnSurfaceEnum : return "NodeOnSurface"; + case NumberNodeToElementConnectivityEnum : return "NumberNodeToElementConnectivity"; case PenaltyOffsetEnum : return "PenaltyOffset"; case PflagEnum : return "Pflag"; Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3946) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3947) @@ -188,4 +188,5 @@ NodeOnIceShelfEnum, NodeOnSurfaceEnum, + NumberNodeToElementConnectivityEnum, PenaltyOffsetEnum, PflagEnum, Index: /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 3946) +++ /issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp (revision 3947) @@ -159,4 +159,5 @@ else if (strcmp(name,"NodeOnIceShelf")==0) return NodeOnIceShelfEnum; else if (strcmp(name,"NodeOnSurface")==0) return NodeOnSurfaceEnum; + else if (strcmp(name,"NumberNodeToElementConnectivity")==0) return NumberNodeToElementConnectivityEnum; else if (strcmp(name,"PenaltyOffset")==0) return PenaltyOffsetEnum; else if (strcmp(name,"Pflag")==0) return PflagEnum; Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 3946) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 3947) @@ -50,41 +50,49 @@ IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); + IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements"); + IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements"); /*2d mesh: */ if (strcmp(iomodel->meshtype,"2d")==0){ - for (i=0;inumberofvertices;i++){ + for (i=0;inumberofelements;i++){ - if(iomodel->my_vertices[i]){ + if(iomodel->my_elements[i]){ - /*Create and add penta element to elements dataset: */ - elements->AddObject(new Sing(i+1,i,iomodel)); + if (*(iomodel->elements_type+2*i+0)==HutterFormulationEnum){ //create only Hutter elements - /*Create and add material property to materials dataset: */ - materials->AddObject(new Matice(i+1,i,iomodel,1)); + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Tria(i+1,i,iomodel)); + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i+1,i,iomodel,3)); + } } - } //for (i=0;inumberofvertices;i++) } //if (strcmp(iomodel->meshtype,"2d")==0) else{ - for (i=0;inumberofvertices;i++){ + for (i=0;inumberofelements;i++){ - if(iomodel->my_vertices[i]){ - if(iomodel->gridonhutter[i]){ - if(!iomodel->gridonsurface[i]){ + if(iomodel->my_elements[i]){ - /*Create and add penta element to elements dataset: */ - elements->AddObject(new Beam(i+1,i,iomodel)); + if (*(iomodel->elements_type+2*i+0)==HutterFormulationEnum){ //create only Hutter elements - /*Create and add material property to materials dataset: */ - materials->AddObject(new Matice(i+1,i,iomodel,2)); + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Penta(i+1,i,iomodel)); - } + + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i+1,i,iomodel,6)); + } } - } //for (i=0;inumberofvertices;i++) + } //for (i=0;inumberofelements;i++) } //if (strcmp(iomodel->meshtype,"2d")==0) @@ -92,20 +100,26 @@ /*Free data: */ xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elements_type); xfree((void**)&iomodel->gridonhutter); xfree((void**)&iomodel->thickness); xfree((void**)&iomodel->surface); xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->gridonsurface); + xfree((void**)&iomodel->gridonbed); + xfree((void**)&iomodel->elementonwater); xfree((void**)&iomodel->uppernodes); xfree((void**)&iomodel->drag_coefficient); xfree((void**)&iomodel->rheology_B); xfree((void**)&iomodel->rheology_n); + xfree((void**)&iomodel->upperelements); + xfree((void**)&iomodel->lowerelements); /*Add new constrant material property to materials, at the end: */ if (strcmp(iomodel->meshtype,"2d")==0){ - materials->AddObject(new Matpar(iomodel->numberofvertices+1,iomodel)); //put it at the end of the materials + materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel)); //put it at the end of the materials } else{ - materials->AddObject(new Matpar(iomodel->numberofvertices2d*(iomodel->numlayers-1)+1,iomodel)); //put it at the end of the materials + materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel)); //put it at the end of the materials } @@ -125,4 +139,6 @@ IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + CreateNumberNodeToElementConnectivity(iomodel); for (i=0;inumberofvertices;i++){ @@ -153,4 +169,6 @@ xfree((void**)&iomodel->gridonicesheet); xfree((void**)&iomodel->gridoniceshelf); + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->numbernodetoelementconnectivity); /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these Index: /issm/trunk/src/c/objects/Elements/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 3946) +++ /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 3947) @@ -414,23 +414,43 @@ int numberofdofspernode; bool onbed; - - + bool onsurface; + int connectivity[2]; + double one0,one1; + + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.deliverp(); + + connectivity[0]=nodes[0]->GetConnectivity(); + connectivity[1]=nodes[1]->GetConnectivity(); + + one0=1/(double)connectivity[0]; + one1=1/(double)connectivity[1]; inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); GetDofList(&doflist[0],&numberofdofspernode); if (onbed){ - Ke_gg[0][0]=1; - Ke_gg[1][1]=1; - Ke_gg[2][0]=-1; - Ke_gg[2][2]=1; - Ke_gg[3][1]=-1; - Ke_gg[3][3]=1; - } - else{ - Ke_gg[2][0]=-1; - Ke_gg[2][2]=1; - Ke_gg[3][1]=-1; - Ke_gg[3][3]=1; + Ke_gg[0][0]=one0; + Ke_gg[1][1]=one0; + Ke_gg[2][0]=-2*one1; + Ke_gg[2][2]=2*one1; + Ke_gg[3][1]=-2*one1; + Ke_gg[3][3]=2*one1; + } + else if (onsurface){ + Ke_gg[2][0]=-one1; + Ke_gg[2][2]=one1; + Ke_gg[3][1]=-one1; + Ke_gg[3][3]=one1; + } + else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity + Ke_gg[2][0]=-2*one1; + Ke_gg[2][2]=2*one1; + Ke_gg[3][1]=-2*one1; + Ke_gg[3][3]=2*one1; } @@ -495,4 +515,6 @@ bool onbed; + bool onsurface; + int connectivity[2]; /*dynamic objects pointed to by hooks: */ @@ -511,4 +533,5 @@ /*recover some inputs: */ inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); /*recover parameters: */ @@ -524,4 +547,7 @@ inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum); + connectivity[0]=nodes[0]->GetConnectivity(); + connectivity[1]=nodes[1]->GetConnectivity(); + //Get all element grid data: GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -551,6 +577,14 @@ GetJacobianDeterminant(&Jdet, &z_list[0],gauss_coord); - for(j=0;jelements); for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab @@ -59,4 +60,5 @@ penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object + ISSMASSERT(iomodel->upperelements); if isnan(iomodel->upperelements[index]){ penta_elements_ids[0]=this->id; //upper penta is the same penta @@ -65,4 +67,5 @@ penta_elements_ids[0]=(int)(iomodel->upperelements[index]); } + ISSMASSERT(iomodel->lowerelements); if isnan(iomodel->lowerelements[index]){ penta_elements_ids[1]=this->id; //lower penta is the same penta @@ -71,5 +74,4 @@ penta_elements_ids[1]=(int)(iomodel->lowerelements[index]); } - this->InitHookNodes(penta_node_ids); this->nodes=NULL; this->InitHookMatice(penta_matice_id);this->matice=NULL; @@ -79,5 +81,4 @@ //intialize inputs, and add as many inputs per element as requested: this->inputs=new Inputs(); - if (iomodel->thickness) { for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1]; @@ -125,5 +126,4 @@ this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs)); } - /*vx,vy and vz: */ if (iomodel->vx) { @@ -178,5 +178,4 @@ this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs)); } - if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); @@ -250,5 +249,4 @@ /*point parameters to real dataset: */ this->parameters=parametersin; - } /*}}}*/ @@ -330,4 +328,6 @@ name==SurfaceEnum || name==BedEnum || + name==SurfaceSlopexEnum || + name==SurfaceSlopeyEnum || name==MeltingRateEnum || name==AccumulationRateEnum || Index: /issm/trunk/src/c/objects/Elements/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 3946) +++ /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 3947) @@ -333,4 +333,16 @@ int doflist[numdofs]; int numberofdofspernode; + int connectivity; + + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.delivers(); + + /*Find connectivity of the node and divide Ke_gg by this connectivity*/ + connectivity=nodes[0]->GetConnectivity(); + Ke_gg[0][0]=1/(double)connectivity; + Ke_gg[1][1]=1/(double)connectivity; GetDofList(&doflist[0],&numberofdofspernode); @@ -374,12 +386,13 @@ double rho_ice,gravity,n,B; double thickness; + int connectivity; /*dynamic objects pointed to by hooks: */ - Node** node=NULL; + Node** nodes=NULL; Matpar* matpar=NULL; Matice* matice=NULL; /*recover objects from hooks: */ - node=(Node**)hnodes.deliverp(); + nodes=(Node**)hnodes.deliverp(); matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); @@ -389,4 +402,7 @@ GetDofList(&doflist[0],&numberofdofspernode); + + //Get connectivity of the node + connectivity=nodes[0]->GetConnectivity(); //compute slope2 @@ -406,9 +422,8 @@ constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); - pe_g[0]=ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0]; - pe_g[1]=vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1]; + pe_g[0]=(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity; + pe_g[1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity; VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES); - } Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3946) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3947) @@ -165,4 +165,6 @@ if (iomodel->gridoniceshelf) this->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,(IssmBool)iomodel->gridoniceshelf[i])); if (iomodel->gridonicesheet) this->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,(IssmBool)iomodel->gridonicesheet[i])); + NumberNodeToElementConnectivityEnum; + if (iomodel->numbernodetoelementconnectivity) this->inputs->AddInput(new IntInput(NumberNodeToElementConnectivityEnum,iomodel->numbernodetoelementconnectivity[i])); } @@ -771,4 +773,14 @@ } /*}}}*/ +/*FUNCTION Node::GetConnectivity {{{2*/ +int Node::GetConnectivity(){ + int connectivity; + + /*recover parameters: */ + inputs->GetParameterValue(&connectivity,NumberNodeToElementConnectivityEnum); + + return connectivity; +} +/*}}}*/ /*FUNCTION Node::GetNumberOfDofs{{{2*/ int Node::GetNumberOfDofs(){ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 3946) +++ /issm/trunk/src/c/objects/Node.h (revision 3947) @@ -70,4 +70,5 @@ int GetDof(int dofindex); void CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s); + int GetConnectivity(); void GetDofList(int* outdoflist,int* pnumberofdofspernode); int GetDofList1(void);