Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3533) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3534) @@ -16,5 +16,7 @@ /*Intermediary*/ - int i,j,k,n; + int i,j; + int vertex_index; + int node_index; /*DataSets: */ @@ -106,10 +108,13 @@ if(iomodel->my_nodes[3*i+j]){ - //Get id of the vertex on which the current node is located - k=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing) - ISSMASSERT(k>=0 && knumberofvertices); + //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(k,i,iomodel)); + nodes->AddObject(new Node(vertex_index,node_index,iomodel)); } Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp (revision 3533) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp (revision 3534) @@ -1,3 +1,3 @@ -/*!\file: CreateParametersPrognostic2.cpp +/*!\file: CreateParametersPrognostic.cpp * \brief driver for creating parameters dataset, for prognostic analysis. */ @@ -9,5 +9,4 @@ #include "../../objects/objects.h" #include "../../shared/shared.h" -#include "../../include/macros.h" #include "../IoModel.h" @@ -17,13 +16,8 @@ DataSet* parameters=NULL; int count; - int i,j; + int i; int dim; - int* part; - - double* vx_g=NULL; - double* vy_g=NULL; - double* h_g=NULL; - double* a_g=NULL; - double* m_g=NULL; + double* vx_g=NULL; + double* vy_g=NULL; /*recover parameters : */ @@ -32,77 +26,89 @@ count=parameters->Size(); - /*Create transformation vector for DG inputs*/ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - part=(int*)xcalloc(iomodel->numberofelements*3,sizeof(int)); - for(i=0;inumberofelements;i++){ - for(j=0;j<3;j++){ - part[3*i+j]=(int)iomodel->elements[3*i+j]-1; //Matlab to C indexing - ISSMASSERT(part[3*i+j]numberofvertices); - } - } - xfree((void**)&iomodel->elements); + /*Get vx and vy: */ + IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); + IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); - /*Get vx: */ - IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); - vx_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(iomodel->vx) for(i=0;i<3*iomodel->numberofelements;i++) vx_g[i]=iomodel->vx[part[i]]/iomodel->yts; + vx_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double)); + vy_g=(double*)xcalloc(iomodel->numberofvertices,sizeof(double)); + + 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,"vx_g",DOUBLEVEC); - param->SetDoubleVec(vx_g,3*iomodel->numberofelements,1); + param->SetDoubleVec(vx_g,iomodel->numberofvertices,1); parameters->AddObject(param); - xfree((void**)&iomodel->vx); - xfree((void**)&vx_g); - - /*Get vy: */ - IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); - vy_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(iomodel->vy) for(i=0;i<3*iomodel->numberofelements;i++) vy_g[i]=iomodel->vy[part[i]]/iomodel->yts; - count++; param= new Param(count,"vy_g",DOUBLEVEC); - param->SetDoubleVec(vy_g,3*iomodel->numberofelements,1); + param->SetDoubleVec(vy_g,iomodel->numberofvertices,1); parameters->AddObject(param); + + xfree((void**)&iomodel->vx); xfree((void**)&iomodel->vy); + xfree((void**)&vx_g); xfree((void**)&vy_g); /*Get thickness: */ IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); - h_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(iomodel->thickness) for(i=0;i<3*iomodel->numberofelements;i++) h_g[i]=iomodel->thickness[part[i]]; - + count++; param= new Param(count,"h_g",DOUBLEVEC); - param->SetDoubleVec(h_g,3*iomodel->numberofelements,1); + 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); - xfree((void**)&h_g); + + /*Get surface: */ + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + + count++; + param= new Param(count,"s_g",DOUBLEVEC); + if(iomodel->surface) param->SetDoubleVec(iomodel->surface,iomodel->numberofvertices,1); + else param->SetDoubleVec(iomodel->surface,0,0); + parameters->AddObject(param); + + /*Free surface: */ + xfree((void**)&iomodel->surface); + + /*Get bed: */ + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + + count++; + param= new Param(count,"b_g",DOUBLEVEC); + if(iomodel->bed) param->SetDoubleVec(iomodel->bed,iomodel->numberofvertices,1); + else param->SetDoubleVec(iomodel->bed,0,0); + parameters->AddObject(param); + + /*Free bed: */ + xfree((void**)&iomodel->bed); + + /*Get melting: */ + IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); + if(iomodel->melting) for(i=0;inumberofvertices;i++)iomodel->melting[i]=iomodel->melting[i]/iomodel->yts; + + count++; + param= new Param(count,"m_g",DOUBLEVEC); + if(iomodel->melting) param->SetDoubleVec(iomodel->melting,iomodel->numberofvertices,1); + else param->SetDoubleVec(iomodel->melting,0,1); + parameters->AddObject(param); + + /*Free melting: */ + xfree((void**)&iomodel->melting); /*Get accumulation: */ IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation"); - a_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(iomodel->accumulation) for(i=0;i<3*iomodel->numberofelements;i++) a_g[i]=iomodel->accumulation[part[i]]/iomodel->yts; - + if(iomodel->accumulation) for(i=0;inumberofvertices;i++)iomodel->accumulation[i]=iomodel->accumulation[i]/iomodel->yts; + count++; param= new Param(count,"a_g",DOUBLEVEC); - param->SetDoubleVec(a_g,3*iomodel->numberofelements,1); + if(iomodel->accumulation) param->SetDoubleVec(iomodel->accumulation,iomodel->numberofvertices,1); + else param->SetDoubleVec(iomodel->accumulation,0,0); parameters->AddObject(param); + + /*Free accumulation: */ xfree((void**)&iomodel->accumulation); - xfree((void**)&a_g); - - /*Get melting: */ - IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); - m_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(iomodel->melting) for(i=0;i<3*iomodel->numberofelements;i++) m_g[i]=iomodel->melting[part[i]]/iomodel->yts; - - count++; - param= new Param(count,"m_g",DOUBLEVEC); - param->SetDoubleVec(m_g,3*iomodel->numberofelements,1); - parameters->AddObject(param); - xfree((void**)&iomodel->melting); - xfree((void**)&m_g); - - /*Free partitioning vector*/ - xfree((void**)&part); /*Assign output pointer: */ Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp =================================================================== --- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 3533) +++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 3534) @@ -22,5 +22,5 @@ /*diverse: */ - int numberofnodes; + int numberofvertices; int i,j,k; int ndof; @@ -32,5 +32,5 @@ /*Some parameters needed: */ - parameters->FindParam(&numberofnodes,"numberofnodes"); + parameters->FindParam(&numberofvertices,"numberofvertices"); /*serialize partition vector: */ @@ -51,10 +51,10 @@ if (ndof!=0){ /*ok, we are dealing with a parameter that needs to be repartitioned, for ndof degrees of freedom: */ - new_parameter=(double*)xcalloc(numberofnodes*ndof,sizeof(double)); + new_parameter=(double*)xcalloc(numberofvertices*ndof,sizeof(double)); param->GetParameterValue(¶meter); for(j=0;jSetDoubleVec(new_parameter,ndof*numberofnodes,ndof); + param->SetDoubleVec(new_parameter,ndof*numberofvertices,ndof); /*Free ressources: */ Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3533) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3534) @@ -580,5 +580,5 @@ for(;;){ - if (node->IsOnSurface())break; + if (node->IsOnSurface()) break; vertexdof=node->GetVertexDof(); @@ -700,5 +700,5 @@ GetDofList(&doflist,&numberofdofspernode); - //initilaize node + //initialize node node=this; @@ -707,5 +707,5 @@ //go through all nodes which sit on top of this node, until we reach the surface, - //and pltg fieldnode in field: + //and plug fieldnode in field: for(;;){ Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3533) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3534) @@ -76,6 +76,16 @@ /*hooks: */ - for(j=0;j<3;j++){ //go recover node ids, needed to initialize the node hook. - tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab + //go recover node ids, needed to initialize the node hook. + if (iomodel->analysis_type==Prognostic2AnalysisEnum()){ + /*Discontinuous Galerkin*/ + tria_node_ids[0]=3*i+1; + tria_node_ids[1]=3*i+2; + tria_node_ids[2]=3*i+3; + } + else{ + /*Continuous Galerkin*/ + for(j=0;j<3;j++){ + tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab + } } tria_matice_id=i+1; //refers to the corresponding ice material object Index: /issm/trunk/src/c/parallel/prognostic2_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic2_core.cpp (revision 3533) +++ /issm/trunk/src/c/parallel/prognostic2_core.cpp (revision 3534) @@ -28,4 +28,5 @@ int numberofdofspernode; int numberofnodes; + int numberofvertices; int dofs[1]={1}; @@ -39,4 +40,5 @@ model->FindParam(&verbose,"verbose"); model->FindParam(&numberofnodes,"numberofnodes"); + model->FindParam(&numberofvertices,"numberofvertices"); model->FindParam(&numberofdofspernode,"numberofdofspernode"); @@ -44,8 +46,10 @@ 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,numberofnodes); - inputs->Add("vy_average",vy_g,1,numberofnodes); + */ + inputs->Add("vx_average",vx_g,1,numberofvertices); + inputs->Add("vy_average",vy_g,1,numberofvertices); _printf_("call computational core:\n"); Index: /issm/trunk/src/m/solutions/jpl/prognostic2.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2.m (revision 3533) +++ /issm/trunk/src/m/solutions/jpl/prognostic2.m (revision 3534) @@ -20,9 +20,9 @@ displaystring(md.verbose,'\n%s',['setup inputs...']); inputs=inputlist; - inputs=add(inputs,'vx',models.p.parameters.vx_g,'doublevec',1,models.p.parameters.numberofnodes); - inputs=add(inputs,'vy',models.p.parameters.vy_g,'doublevec',1,models.p.parameters.numberofnodes); - inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes); - inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes); - inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.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,'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); inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double'); Index: /issm/trunk/src/m/solutions/jpl/prognostic2_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 3533) +++ /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 3534) @@ -14,8 +14,11 @@ vx_g=get(inputs,'vx',1); vy_g=get(inputs,'vy',1); - 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.numberofnodes); - inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.numberofnodes); + + %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:']);