Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp (revision 3582) @@ -11,36 +11,9 @@ void CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){ - - int i; - int count=0; + DataSet* constraints = NULL; - + /*Create constraints: */ constraints = new DataSet(ConstraintsEnum); - - /*Fetch data: */ - IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,"spcthickness"); - - count=1; //matlab indexing - /*Create spcs from x,y,z, as well as the spc values on those spcs: */ - for (i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - - if ((int)iomodel->spcthickness[2*i]){ - - constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcthickness+2*i+1)));//we enforce first translation degree of freedom, for temperature - count++; - } - } - } - - /*Free data: */ - xfree((void**)&iomodel->spcthickness); - - /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these - * datasets, it will not be redone: */ - constraints->Presort(); - - cleanup_and_return: /*Assign output pointer: */ Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (revision 3582) @@ -8,18 +8,21 @@ #include "../../objects/objects.h" #include "../../shared/shared.h" +#include "../../include/macros.h" #include "../../include/typedefs.h" #include "../../MeshPartitionx/MeshPartitionx.h" #include "../IoModel.h" -void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ /*Intermediary*/ - int i,j,k,n; + int i,j; + int vertex_index; + int node_index; /*DataSets: */ - DataSet* elements = NULL; - DataSet* nodes = NULL; - DataSet* vertices = NULL; - DataSet* materials = NULL; + DataSet* elements = NULL; + DataSet* nodes = NULL; + DataSet* vertices = NULL; + DataSet* materials = NULL; /*First create the elements, nodes and material properties: */ @@ -32,4 +35,5 @@ Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ /*2d mesh: */ if (strcmp(iomodel->meshtype,"2d")==0){ @@ -54,7 +58,6 @@ } }//for (i=0;ielements); xfree((void**)&iomodel->thickness); xfree((void**)&iomodel->surface); @@ -65,45 +68,11 @@ } else{ // if (strcmp(meshtype,"2d")==0) - - /*Fetch data needed: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); - IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); - IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); - IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); - IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); - IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); - IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); - - for (i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - /*Create and add penta element to elements dataset: */ - elements->AddObject(new Penta(i,iomodel)); - - /*Create and add material property to materials dataset: */ - materials->AddObject(new Matice(i,iomodel,6)); - } - }//for (i=0;ielements); - xfree((void**)&iomodel->thickness); - xfree((void**)&iomodel->surface); - xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->elementoniceshelf); - xfree((void**)&iomodel->elementonbed); - xfree((void**)&iomodel->elementonsurface); - xfree((void**)&iomodel->elementonwater); - + ISSMERROR("not implemented yet"); } //if (strcmp(meshtype,"2d")==0) - /*Add new constrant material property to materials, at the end: */ + /*Add new constrant material property tgo materials, at the end: */ materials->AddObject(new Matpar(iomodel)); - /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ - IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); - IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); - } + /*Create nodes and vertices: */ IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); @@ -113,7 +82,13 @@ IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); + IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + if (strcmp(iomodel->meshtype,"3d")==0){ + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); + } + /*Build Vertices dataset*/ for (i=0;inumberofvertices;i++){ @@ -124,7 +99,24 @@ vertices->AddObject(new Vertex(i,iomodel)); - /*Add node to nodes dataset: */ - nodes->AddObject(new Node(i,iomodel)); + } + } + /*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/ + for (i=0;inumberofelements;i++){ + for (j=0;j<3;j++){ + + if(iomodel->my_nodes[3*i+j]){ + + //Get index of the vertex on which the current node is located + vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing) + ISSMASSERT(vertex_index>=0 && vertex_indexnumberofvertices); + + //Compute Node index (id-1) + node_index=3*i+j; + + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(vertex_index,node_index,iomodel)); + + } } } @@ -139,4 +131,5 @@ xfree((void**)&iomodel->gridonbed); xfree((void**)&iomodel->gridonsurface); + xfree((void**)&iomodel->gridonhutter); xfree((void**)&iomodel->uppernodes); xfree((void**)&iomodel->gridonicesheet); @@ -146,6 +139,6 @@ * datasets, it will not be redone: */ elements->Presort(); + nodes->Presort(); vertices->Presort(); - nodes->Presort(); materials->Presort(); @@ -157,4 +150,3 @@ *pvertices=vertices; *pmaterials=materials; - } Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp (revision 3582) @@ -8,20 +8,108 @@ #include "../../shared/shared.h" #include "../../include/macros.h" +#include "../../include/typedefs.h" #include "../IoModel.h" - void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){ - DataSet* loads = NULL; + /*Intermediary*/ + int i,j; + int i1,i2; + int pos1,pos2; + double e1,e2; + + /*Output*/ + DataSet* loads=NULL; + + /*numericalflux intermediary data: */ + char numericalflux_type[NUMERICALFLUXSTRING]; + int numericalflux_id; + int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES]; + int numericalflux_elem_id; + double numericalflux_h[MAX_NUMERICALFLUX_NODES]; /*Create loads: */ loads = new DataSet(LoadsEnum); - - + + /*Get edges and elements*/ + IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges"); + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + + /*First load data:*/ + for (i=0;inumberofedges;i++){ + + /*Get left and right elements*/ + e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] + e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] + + /*Now, if this element is not in the partition, pass: */ + if(!iomodel->my_elements[(int)e1]) continue; + + /*Create load*/ + numericalflux_id=i+1; //Matlab indexing + numericalflux_elem_id=(int)e1+1;//id is in matlab index + + /*1: Get vertices ids*/ + i1=(int)iomodel->edges[4*i+0]; + i2=(int)iomodel->edges[4*i+1]; + + if (!isnan(e2)){ + strcpy(numericalflux_type,"internal"); + + /*Now, we must get the nodes of the 4 nodes located on the edge*/ + + /*2: Get the column where these ids are located in the index*/ + pos1=pos2=UNDEF; + for(j=0;j<3;j++){ + if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1; + if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1; + } + ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF); + + /*3: We have the id of the elements and the position of the vertices in the index + * we can compute their dofs!*/ + numericalflux_node_ids[0]=3*(int)e1+pos1; //ex: 1 2 3 + numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1 + numericalflux_node_ids[2]=3*(int)e2+pos2; //ex: 1 2 3 + numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2 + + numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1] -1]; + numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1]; + numericalflux_h[2]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[2]-1]-1]; + numericalflux_h[3]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[3]-1]-1]; + } + else{ + strcpy(numericalflux_type,"boundary"); + + /*2: Get the column where these ids are located in the index*/ + pos1==UNDEF; + for(j=0;j<3;j++){ + if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1; + } + ISSMASSERT(pos1!=UNDEF); + + /*3: We have the id of the elements and the position of the vertices in the index + * we can compute their dofs!*/ + numericalflux_node_ids[0]=3*(int)e1+pos1; + numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; + + numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1]-1]; + numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1]; + numericalflux_h[2]=UNDEF; + numericalflux_h[3]=UNDEF; + } + + loads->AddObject(new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id,numericalflux_h)); + } + + /*Free data: */ + xfree((void**)&iomodel->edges); + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->thickness); + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these * datasets, it will not be redone: */ loads->Presort(); - - cleanup_and_return: /*Assign output pointer: */ @@ -29,4 +117,2 @@ } - - Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp (revision 3582) @@ -1,3 +1,4 @@ - /* \brief driver for creating parameters dataset, for prognostic analysis. +/*!\file: CreateParametersPrognostic.cpp + * \brief driver for creating parameters dataset, for prognostic analysis. */ @@ -15,5 +16,7 @@ int count; int i; - double* u_g=NULL; + int dim; + double* vx_g=NULL; + double* vy_g=NULL; /*recover parameters : */ @@ -25,22 +28,49 @@ IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); - IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz"); - u_g=(double*)xcalloc(iomodel->numberofvertices*3,sizeof(double)); + vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double)); + vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double)); - if(iomodel->vx) for(i=0;inumberofvertices;i++)u_g[3*i+0]=iomodel->vx[i]/iomodel->yts; - if(iomodel->vy) for(i=0;inumberofvertices;i++)u_g[3*i+1]=iomodel->vy[i]/iomodel->yts; - if(iomodel->vz) for(i=0;inumberofvertices;i++)u_g[3*i+2]=iomodel->vz[i]/iomodel->yts; + if(iomodel->vx) for(i=0;inumberofvertices;i++)vx_g[i]=iomodel->vx[i]/iomodel->yts; + if(iomodel->vy) for(i=0;inumberofvertices;i++)vy_g[i]=iomodel->vy[i]/iomodel->yts; count++; - param= new Param(count,"u_g",DOUBLEVEC); - param->SetDoubleVec(u_g,3*iomodel->numberofvertices,3); + param= new Param(count,"vx_g",DOUBLEVEC); + param->SetDoubleVec(vx_g,iomodel->numberofvertices,1); parameters->AddObject(param); - + count++; + param= new Param(count,"vy_g",DOUBLEVEC); + param->SetDoubleVec(vy_g,iomodel->numberofvertices,1); + parameters->AddObject(param); xfree((void**)&iomodel->vx); xfree((void**)&iomodel->vy); - xfree((void**)&iomodel->vz); - xfree((void**)&u_g); + xfree((void**)&vx_g); + xfree((void**)&vy_g); + + /*Get thickness: */ + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + + count++; + param= new Param(count,"h_g",DOUBLEVEC); + if(iomodel->thickness) param->SetDoubleVec(iomodel->thickness,iomodel->numberofvertices,1); + else param->SetDoubleVec(iomodel->thickness,0,0); + parameters->AddObject(param); + + /*Free thickness: */ + xfree((void**)&iomodel->thickness); + + /*Get dhdt: */ + IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt"); + if(iomodel->dhdt) for(i=0;inumberofvertices;i++)iomodel->dhdt[i]=iomodel->dhdt[i]/iomodel->yts; + + count++; + param= new Param(count,"dhdt_g",DOUBLEVEC); + if(iomodel->dhdt) param->SetDoubleVec(iomodel->dhdt,iomodel->numberofvertices,1); + else param->SetDoubleVec(iomodel->dhdt,0,1); + parameters->AddObject(param); + + /*Free dhdt: */ + xfree((void**)&iomodel->dhdt); /*Get melting: */ Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3582) @@ -245,5 +245,9 @@ count++; param= new Param(count,"numberofnodes",DOUBLE); - if (iomodel->analysis_type==Prognostic2AnalysisEnum) param->SetDouble(3*iomodel->numberofelements); + if ( + iomodel->analysis_type==Prognostic2AnalysisEnum || + iomodel->analysis_type==BalancedthicknessAnalysisEnum + ) + param->SetDouble(3*iomodel->numberofelements); else param->SetDouble(iomodel->numberofvertices); parameters->AddObject(param); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3582) @@ -171,4 +171,5 @@ /*!basal: */ iomodel->accumulation=NULL; + iomodel->dhdt=NULL; /*parameter output: */ @@ -254,4 +255,5 @@ xfree((void**)&iomodel->melting); xfree((void**)&iomodel->accumulation); + xfree((void**)&iomodel->dhdt); xfree((void**)&iomodel->B); xfree((void**)&iomodel->n); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3582) @@ -172,4 +172,5 @@ double* melting; double* accumulation; + double* dhdt; /*parameter output: */ Index: /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp (revision 3581) +++ /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp (revision 3582) @@ -22,5 +22,5 @@ void Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ - if (iomodel->analysis_type==Prognostic2AnalysisEnum) + if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum) DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle); else Index: /issm/trunk/src/c/objects/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3581) +++ /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3582) @@ -222,6 +222,14 @@ found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); if(!found)ISSMERROR(" could not find vy_average in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + if (analysis_type==Prognostic2AnalysisEnum){ + if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + /*No transient term is involved*/ + dt=1; + } + else{ + ISSMERROR("analysis_type %i not supported yet"); + } /* Get node coordinates, dof list and normal vector: */ @@ -330,6 +338,14 @@ found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); if(!found)ISSMERROR(" could not find vy_average in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + if (analysis_type==Prognostic2AnalysisEnum){ + if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + /*No transient term is involved*/ + dt=1; + } + else{ + ISSMERROR("analysis_type %i not supported yet"); + } /* Get node coordinates, dof list and normal vector: */ @@ -464,6 +480,14 @@ found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); if(!found)ISSMERROR(" could not find vy_average in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + if (analysis_type==Prognostic2AnalysisEnum){ + if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + /*No transient term is involved*/ + dt=1; + } + else{ + ISSMERROR("analysis_type %i not supported yet"); + } /* Get node coordinates, dof list and normal vector: */ Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3581) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3582) @@ -585,12 +585,4 @@ CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type); } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - - CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - - CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type); - } else if (analysis_type==ThermalAnalysisEnum){ @@ -604,46 +596,4 @@ ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet"); } - -} -/*}}}*/ -/*FUNCTION CreateKMatrixBalancedthickness {{{1*/ - -void Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(this->properties.onwater)return; - - /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; - return; - -} -/*}}}*/ -/*FUNCTION CreateKMatrixBalancedvelocities {{{1*/ - -void Penta::CreateKMatrixBalancedvelocities(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(this->properties.onwater)return; - - /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; - return; } @@ -1711,12 +1661,4 @@ CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - - CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - - CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type); - } else if (analysis_type==ThermalAnalysisEnum){ @@ -1731,44 +1673,4 @@ } -} -/*}}}*/ -/*FUNCTION CreatePVectorBalancedthickness {{{1*/ - -void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(this->properties.onwater)return; - - /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); - delete tria; - return; -} -/*}}}*/ -/*FUNCTION CreatePVectorBalancedvelocities {{{1*/ - -void Penta::CreatePVectorBalancedvelocities( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(this->properties.onwater)return; - - /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); - delete tria; - return; } /*}}}*/ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3581) +++ /issm/trunk/src/c/objects/Penta.h (revision 3582) @@ -110,8 +110,4 @@ void CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); void CreatePVectorPrognostic( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); void CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3581) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3582) @@ -77,5 +77,5 @@ /*hooks: */ //go recover node ids, needed to initialize the node hook. - if (iomodel->analysis_type==Prognostic2AnalysisEnum){ + if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==BalancedthicknessAnalysisEnum){ /*Discontinuous Galerkin*/ tria_node_ids[0]=3*i+1; @@ -633,4 +633,121 @@ /* matrices: */ + double B[2][numgrids]; + double Bprime[2][numgrids]; + double DL[2][2]={0.0}; + double DLprime[2][2]={0.0}; + double DL_scalar; + double Ke_gg[numdof][numdof]={0.0}; + double Ke_gg2[numdof][numdof]={0.0}; + double Jdettria; + + /*input parameters for structural analysis (diagnostic): */ + double vx_list[numgrids]={0.0}; + double vy_list[numgrids]={0.0}; + double vx,vy; + int dofs[1]={0}; + int found; + + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + Matpar* matpar=NULL; + Matice* matice=NULL; + Numpar* numpar=NULL; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover objects from hooks: */ + nodes=(Node**) hnodes.deliverp(); + matpar=(Matpar*)hmatpar.delivers(); + matice=(Matice*)hmatice.delivers(); + numpar=(Numpar*)hnumpar.delivers(); + + /*recover extra inputs from users, at current convergence iteration: */ + found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vx_average in inputs!"); + found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vy_average in inputs!"); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igproperties.s[i]; + + /*Get normal vector to the surface*/ + nx=(vx_list[0]+vx_list[1]+vx_list[2])/3; + ny=(vy_list[0]+vy_list[1]+vy_list[2])/3; + if(nx==0 && ny==0){ + SurfaceNormal(&surface_normal[0],xyz_list); + nx=surface_normal[0]; + ny=surface_normal[1]; + } + if(nx==0 && ny==0){ + nx=0; + ny=1; + } + norm=pow( pow(nx,2)+pow(ny,2) , (double).5); + nx=nx/norm; + ny=ny/norm; + //Create Artificial diffusivity once for all if requested if(numpar->artdiff){ @@ -700,6 +837,6 @@ v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); - K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); - K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! + K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); } @@ -734,27 +871,16 @@ DL_scalar=gauss_weight*Jdettria; - //Create DL and DLprime matrix - DL[0][0]=DL_scalar*dvxdx; - DL[1][1]=DL_scalar*dvydy; - - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; + DLprime[0][0]=DL_scalar*nx; + DLprime[1][1]=DL_scalar*ny; //Do the triple product tL*D*L. - //Ke_gg_thickness=B'*DLprime*Bprime; - - TripleMultiply( &B[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke_gg_thickness1[0][0],0); - + //Ke_gg_velocities=B'*DLprime*Bprime; TripleMultiply( &B[0][0],2,numdof,1, &DLprime[0][0],2,2,0, &Bprime[0][0],2,numdof,0, - &Ke_gg_thickness2[0][0],0); + &Ke_gg_velocities2[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iartdiff){ @@ -822,6 +948,6 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/ -void Tria::CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +/*FUNCTION Tria::CreateKMatrixDiagnosticHoriz {{{1*/ +void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -830,6 +956,5 @@ /* node data: */ const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; + const int numdof=2*numgrids; double xyz_list[numgrids][3]; int doflist[numdof]; @@ -845,4 +970,564 @@ double gauss_l1l2l3[3]; + /* material data: */ + double viscosity; //viscosity + double newviscosity; //viscosity + double oldviscosity; //viscosity + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/ + + /* matrices: */ + double B[3][numdof]; + double Bprime[3][numdof]; + double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. + double D_scalar; + + /* local element matrices: */ + double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. + + double Jdet; + + /*input parameters for structural analysis (diagnostic): */ + double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; + double oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; + double thickness; + int dofs[2]={0,1}; + + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + Matpar* matpar=NULL; + Matice* matice=NULL; + Numpar* numpar=NULL; + + ParameterInputs* inputs=NULL; + + /*First, if we are on water, return empty matrix: */ + if(this->properties.onwater) return; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover objects from hooks: */ + nodes =(Node**) hnodes.deliverp(); + matpar=(Matpar*)hmatpar.delivers(); + matice=(Matice*)hmatice.delivers(); + numpar=(Numpar*)hnumpar.delivers(); + + /*recover extra inputs from users, at current convergence iteration: */ + inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iproperties.h[0],gauss_l1l2l3); + + /*Get strain rate from velocity: */ + GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); + GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get viscosity: */ + matice->GetViscosity2d(&viscosity, &epsilon[0]); + matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=newviscosity*thickness*gauss_weight*Jdet; + + for (i=0;i<3;i++){ + D[i][i]=D_scalar; + } + + /*Get B and Bprime matrices: */ + GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],3,numdof,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iproperties.shelf){ + CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type); + } + +cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/ +void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ + + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int numdof=2*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* matrices: */ + double L[2][numdof]; + double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + double DL_scalar; + + /* local element matrices: */ + double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag + + double Jdet; + + /*slope: */ + double slope[2]={0.0,0.0}; + double slope_magnitude; + + /*input parameters for structural analysis (diagnostic): */ + double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; + int dofs[2]={0,1}; + + /*friction: */ + double alpha2_list[numgrids]={0.0,0.0,0.0}; + double alpha2; + + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + Matpar* matpar=NULL; + Matice* matice=NULL; + Numpar* numpar=NULL; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.deliverp(); + matpar=(Matpar*)hmatpar.delivers(); + matice=(Matice*)hmatice.delivers(); + numpar=(Numpar*)hnumpar.delivers(); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iproperties.shelf){ + /*no friction, do nothing*/ + return; + } + + if (this->properties.friction_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); + + /*recover extra inputs from users, at current convergence iteration: */ + inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + + /*Build alpha2_list used by drag stiffness matrix*/ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=matpar->GetG(); + friction->rho_ice=matpar->GetRhoIce(); + friction->rho_water=matpar->GetRhoWater(); + friction->K=&this->properties.k[0]; + friction->bed=&this->properties.b[0]; + friction->thickness=&this->properties.h[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=this->properties.p; + friction->q=this->properties.q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); + } + #endif + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, + //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ + GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_l1l2l3); + slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); + + if (slope_magnitude>MAXSLOPE){ + alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); + alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); + alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); + } + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get L matrix: */ + GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + + /*Now, take care of the basal friction if there is any: */ + GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); + + DL_scalar=alpha2*gauss_weight*Jdet; + for (i=0;i<2;i++){ + DL[i][i]=DL_scalar; + } + + /* Do the triple producte tL*D*L: */ + TripleMultiply( &L[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &L[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + for( i=0; iGetLatentHeat(); + heatcapacity=matpar->GetHeatCapacity(); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights: */ + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igRecover("dt",&dt); + if(!found)ISSMERROR(" could not find dt in inputs!"); + /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - /*Modify z so that it reflects the surface*/ - for(i=0;iproperties.s[i]; - - /*Get normal vector to the surface*/ - nx=(vx_list[0]+vx_list[1]+vx_list[2])/3; - ny=(vy_list[0]+vy_list[1]+vy_list[2])/3; - if(nx==0 && ny==0){ - SurfaceNormal(&surface_normal[0],xyz_list); - nx=surface_normal[0]; - ny=surface_normal[1]; - } - if(nx==0 && ny==0){ - nx=0; - ny=1; - } - norm=pow( pow(nx,2)+pow(ny,2) , (double).5); - nx=nx/norm; - ny=ny/norm; - //Create Artificial diffusivity once for all if requested - if(numpar->artdiff){ + if(numpar->artdiff==1){ //Get the Jacobian determinant gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; @@ -933,6 +1601,6 @@ v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); - K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! - K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); + K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); } @@ -942,4 +1610,5 @@ /* Start looping on the number of gaussian points: */ for (ig=0; igartdiff){ + for( i=0; iartdiff==1){ /* Compute artificial diffusivity */ @@ -1044,792 +1736,4 @@ } /*}}}*/ -/*FUNCTION Tria::CreateKMatrixDiagnosticHoriz {{{1*/ -void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=3; - const int numdof=2*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - /* material data: */ - double viscosity; //viscosity - double newviscosity; //viscosity - double oldviscosity; //viscosity - - /* strain rate: */ - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/ - - /* matrices: */ - double B[3][numdof]; - double Bprime[3][numdof]; - double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. - double D_scalar; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - - double Jdet; - - /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - double oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - double thickness; - int dofs[2]={0,1}; - - /*dynamic objects pointed to by hooks: */ - Node** nodes=NULL; - Matpar* matpar=NULL; - Matice* matice=NULL; - Numpar* numpar=NULL; - - ParameterInputs* inputs=NULL; - - /*First, if we are on water, return empty matrix: */ - if(this->properties.onwater) return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover objects from hooks: */ - nodes =(Node**) hnodes.deliverp(); - matpar=(Matpar*)hmatpar.delivers(); - matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;iproperties.h[0],gauss_l1l2l3); - - /*Get strain rate from velocity: */ - GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); - GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); - - /*Get viscosity: */ - matice->GetViscosity2d(&viscosity, &epsilon[0]); - matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ - newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*thickness*gauss_weight*Jdet; - - for (i=0;i<3;i++){ - D[i][i]=D_scalar; - } - - /*Get B and Bprime matrices: */ - GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iproperties.shelf){ - CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type); - } - -cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/ -void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ - - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=3; - const int numdof=2*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - /* matrices: */ - double L[2][numdof]; - double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag - double DL_scalar; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag - - double Jdet; - - /*slope: */ - double slope[2]={0.0,0.0}; - double slope_magnitude; - - /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - int dofs[2]={0,1}; - - /*friction: */ - double alpha2_list[numgrids]={0.0,0.0,0.0}; - double alpha2; - - double MAXSLOPE=.06; // 6 % - double MOUNTAINKEXPONENT=10; - - /*dynamic objects pointed to by hooks: */ - Node** nodes=NULL; - Matpar* matpar=NULL; - Matice* matice=NULL; - Numpar* numpar=NULL; - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover objects from hooks: */ - nodes=(Node**)hnodes.deliverp(); - matpar=(Matpar*)hmatpar.delivers(); - matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;iproperties.shelf){ - /*no friction, do nothing*/ - return; - } - - if (this->properties.friction_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - - /*Build alpha2_list used by drag stiffness matrix*/ - Friction* friction=NewFriction(); - - /*Initialize all fields: */ - friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); - strcpy(friction->element_type,"2d"); - - friction->gravity=matpar->GetG(); - friction->rho_ice=matpar->GetRhoIce(); - friction->rho_water=matpar->GetRhoWater(); - friction->K=&this->properties.k[0]; - friction->bed=&this->properties.b[0]; - friction->thickness=&this->properties.h[0]; - friction->velocities=&vxvy_list[0][0]; - friction->p=this->properties.p; - friction->q=this->properties.q; - - /*Compute alpha2_list: */ - FrictionGetAlpha2(&alpha2_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); - } - #endif - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, - //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ - GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_l1l2l3); - slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - - if (slope_magnitude>MAXSLOPE){ - alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); - } - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - - /*Now, take care of the basal friction if there is any: */ - GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); - - DL_scalar=alpha2*gauss_weight*Jdet; - for (i=0;i<2;i++){ - DL[i][i]=DL_scalar; - } - - /* Do the triple producte tL*D*L: */ - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; iGetLatentHeat(); - heatcapacity=matpar->GetHeatCapacity(); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Get gaussian points and weights: */ - GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igRecover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity_average in inputs!"); - - for(i=0;iRecover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - //Create Artificial diffusivity once for all if requested - if(numpar->artdiff==1){ - //Get the Jacobian determinant - gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); - - //Build K matrix (artificial diffusivity matrix) - v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); - v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); - - K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); - K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igartdiff==1){ - - /* Compute artificial diffusivity */ - KDL[0][0]=DL_scalar*K[0][0]; - KDL[1][1]=DL_scalar*K[1][1]; - - TripleMultiply( &Bprime[0][0],2,numdof,1, - &KDL[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add artificial diffusivity matrix */ - for( i=0; iRecover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); if(!found)ISSMERROR(" could not find melting in inputs!"); + found=inputs->Recover("dhdt",&dhdt_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find dhdt in inputs!"); /* Get node coordinates and dof list: */ @@ -2301,7 +2207,8 @@ GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); + GetParameterValue(&dhdt_g, &dhdt_list[0],gauss_l1l2l3); /* Add value into pe_g: */ - for( i=0; iFindParam(&u_g_serial,NULL,NULL,"u_g",BalancedthicknessAnalysisEnum); - model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum); - model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum); + model->FindParam(&vx_g,NULL,NULL,"vx_g",BalancedthicknessAnalysisEnum); + model->FindParam(&vy_g,NULL,NULL,"vy_g",BalancedthicknessAnalysisEnum); + model->FindParam(&m_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum); + model->FindParam(&a_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum); + model->FindParam(&h_g,NULL,NULL,"h_g",BalancedthicknessAnalysisEnum); + model->FindParam(&dhdt_g,NULL,NULL,"dhdt_g",BalancedthicknessAnalysisEnum); model->FindParam(&numberofnodes,"numberofnodes"); inputs=new ParameterInputs; - inputs->Add("velocity",u_g_serial,3,numberofnodes); - inputs->Add("melting",melting_g,1,numberofnodes); - inputs->Add("accumulation",accumulation_g,1,numberofnodes); + inputs->Add("vx",vx_g,1,numberofnodes); + inputs->Add("vy",vy_g,1,numberofnodes); + inputs->Add("melting",m_g,1,numberofnodes); + inputs->Add("accumulation",a_g,1,numberofnodes); + inputs->Add("thickness",h_g,1,numberofnodes); + inputs->Add("dhdt",dhdt_g,1,numberofnodes); _printf_("initialize results:\n"); @@ -138,6 +146,13 @@ /*Free ressources:*/ + xfree((void**)&vx_g); + xfree((void**)&vy_g); + xfree((void**)&m_g); + xfree((void**)&a_g); + xfree((void**)&dhdt_g); + xfree((void**)&h_g); delete processedresults; delete results; + delete model; delete inputs; Index: /issm/trunk/src/c/parallel/balancedthickness_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/balancedthickness_core.cpp (revision 3581) +++ /issm/trunk/src/c/parallel/balancedthickness_core.cpp (revision 3582) @@ -19,5 +19,6 @@ /*intermediary: */ - Vec u_g=NULL; + Vec vx_g=NULL; + Vec vy_g=NULL; /*solutions: */ @@ -28,5 +29,6 @@ int numberofdofspernode; int numberofnodes; - int dofs[2]={1,1}; + int numberofvertices; + int dofs[1]={1}; /*fem balancedthickness model: */ @@ -34,5 +36,4 @@ - /*recover fem model: */ fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum); @@ -40,16 +41,25 @@ model->FindParam(&verbose,"verbose"); model->FindParam(&numberofnodes,"numberofnodes"); + model->FindParam(&numberofvertices,"numberofvertices"); model->FindParam(&numberofdofspernode,"numberofdofspernode"); _printf_("depth averaging velocity...\n"); - u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity - FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); - inputs->Add("velocity_average",u_g,2,numberofnodes); + vx_g=inputs->Get("vx",&dofs[0],1); + vy_g=inputs->Get("vy",&dofs[0],1); + /* NOT WORKING YET.... + FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx"); + FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy"); + */ + inputs->Add("vx_average",vx_g,1,numberofvertices); + inputs->Add("vy_average",vy_g,1,numberofvertices); _printf_("call computational core:\n"); diagnostic_core_linear(&h_g,fem_p,inputs,BalancedthicknessAnalysisEnum,NoneAnalysisEnum); - _printf_("extrude computed thickness on all layers:\n"); - FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); + _printf_("Averaging over vertices:\n"); + FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters); + +// _printf_("extrude computed thickness on all layers:\n"); +// FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); /*Plug results into output dataset: */ @@ -58,5 +68,6 @@ /*Free ressources:*/ - VecFree(&u_g); + VecFree(&vx_g); + VecFree(&vy_g); VecFree(&h_g); } Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 3581) +++ /issm/trunk/src/m/classes/@model/model.m (revision 3582) @@ -161,4 +161,5 @@ md.vel_obs_raw=NaN; md.accumulation=NaN; + md.dhdt=NaN; md.geothermalflux=NaN; md.observed_temperature=NaN; Index: /issm/trunk/src/m/classes/public/display/displayobservations.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displayobservations.m (revision 3581) +++ /issm/trunk/src/m/classes/public/display/displayobservations.m (revision 3582) @@ -18,4 +18,5 @@ fielddisplay(md,'vel_obs_raw','raw observed magnitude [m/a]'); fielddisplay(md,'accumulation','surface accumulation rate [m/a]'); +fielddisplay(md,'dhdt','surface dhdt rate [m/a]'); fielddisplay(md,'observed_temperature','observed temperature [K]'); fielddisplay(md,'geothermalflux','geothermal heat flux [W/m^2]'); Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 3581) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 3582) @@ -340,12 +340,8 @@ %VELOCITIES MELTING AND ACCUMULATION - fields={'vx','vy','accumulation','melting'}; - checksize(md,fields,[md.numberofgrids 1]); - checknan(md,fields); - - %SPC - if any(md.spcthickness(find(md.gridonboundary))~=1), - error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']); - end + fields={'vx','vy','accumulation','melting','dhdt'}; + checksize(md,fields,[md.numberofgrids 1]); + checknan(md,fields); + end Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 3581) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 3582) @@ -87,4 +87,5 @@ WriteData(fid,md.melting,'Mat','melting'); WriteData(fid,md.accumulation,'Mat','accumulation'); +WriteData(fid,md.dhdt,'Mat','dhdt'); %Get materials Index: /issm/trunk/src/m/solutions/jpl/balancedthickness.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedthickness.m (revision 3581) +++ /issm/trunk/src/m/solutions/jpl/balancedthickness.m (revision 3582) @@ -12,5 +12,5 @@ displaystring(md.verbose,'%s',['reading balancedthickness model data']); - md.analysis_type=BalancedthicknessAnalysisEnum; models.bt=CreateFemModel(md); + md.analysis_type=BalancedthicknessAnalysisEnum; models.p=CreateFemModel(md); % figure out number of dof: just for information purposes. @@ -20,7 +20,10 @@ displaystring(md.verbose,'\n%s',['setup inputs...']); inputs=inputlist; - inputs=add(inputs,'velocity',models.bt.parameters.u_g,'doublevec',3,models.bt.parameters.numberofnodes); - inputs=add(inputs,'melting',models.bt.parameters.m_g,'doublevec',1,models.bt.parameters.numberofnodes); - inputs=add(inputs,'accumulation',models.bt.parameters.a_g,'doublevec',1,models.bt.parameters.numberofnodes); + inputs=add(inputs,'vx',models.p.parameters.vx_g,'doublevec',1,models.p.parameters.numberofvertices); + inputs=add(inputs,'vy',models.p.parameters.vy_g,'doublevec',1,models.p.parameters.numberofvertices); + inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofvertices); + inputs=add(inputs,'dhdt',models.p.parameters.dhdt_g,'doublevec',1,models.p.parameters.numberofvertices); + inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofvertices); + inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofvertices); displaystring(md.verbose,'\n%s',['call computational core:']); Index: /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m (revision 3581) +++ /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m (revision 3582) @@ -6,5 +6,5 @@ %get FE model - m=models.bt; + m=models.p; results.time=0; results.step=1; @@ -12,13 +12,22 @@ displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']); %Take only the first two dofs of m.parameters.u_g - u_g=get(inputs,'velocity',[1 1 0 0]); - u_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,u_g,'velocity'); - inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes); + vx_g=get(inputs,'vx',1); + vy_g=get(inputs,'vy',1); + + %NOT WORKING YET!!!!! + %vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx'); + %vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy'); + + inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.numberofvertices); + inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.numberofvertices); displaystring(m.parameters.verbose,'\n%s',['call computational core:']); results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type); + displaystring(m.parameters.verbose,'\n%s',['averaging over vertices']); + results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g); + displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']); - results.h_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g,'thickness',0); + %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness'); end %end function Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m (revision 3581) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m (revision 3582) @@ -13,4 +13,5 @@ %system matrices [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + save A K_gg p_g [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/m/utils/BC/SetIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetIceSheetBC.m (revision 3581) +++ /issm/trunk/src/m/utils/BC/SetIceSheetBC.m (revision 3582) @@ -36,4 +36,8 @@ disp(' no melting specified: values set as zero'); end +if isnan(md.dhdt), + md.dhdt=zeros(md.numberofgrids,1); + disp(' no dhdt specified: values set as zero'); +end displaystring(md.verbose,'%s',[' boundary conditions for prognostic model initialization']); Index: /issm/trunk/src/m/utils/BC/SetIceShelfBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetIceShelfBC.m (revision 3581) +++ /issm/trunk/src/m/utils/BC/SetIceShelfBC.m (revision 3582) @@ -60,4 +60,8 @@ disp(' no melting specified: values set as zero'); end +if isnan(md.dhdt), + md.dhdt=zeros(md.numberofgrids,1); + disp(' no dhdt specified: values set as zero'); +end displaystring(md.verbose,'%s',[' boundary conditions for prognostic model initialization']); Index: /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m =================================================================== --- /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 3581) +++ /issm/trunk/src/m/utils/BC/SetMarineIceSheetBC.m (revision 3582) @@ -71,4 +71,8 @@ disp(' no melting specified: values set as zero'); end +if isnan(md.dhdt), + md.dhdt=zeros(md.numberofgrids,1); + disp(' no dhdt specified: values set as zero'); +end displaystring(md.verbose,'%s',[' boundary conditions for prognostic model initialization']);