Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24087) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24088) @@ -2021,4 +2021,14 @@ } _error_("Node provided not found among element nodes"); +} +/*}}}*/ +int Tria::GetVertexIndex(Vertex* vertex){/*{{{*/ + + _assert_(vertices); + for(int i=0;i(numnodes); + tria_node_ids[0]= index + 1; + break; case P1Enum: numnodes = 3; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24087) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24088) @@ -85,4 +85,5 @@ void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level); int GetNodeIndex(Node* node); + int GetVertexIndex(Vertex* vertex); int GetNumberOfNodes(void); int GetNumberOfNodes(int enum_type); Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24087) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24088) @@ -17,6 +17,4 @@ /*Load macros*/ #define NUMVERTICES 2 -#define NUMNODES_INTERNAL 4 -#define NUMNODES_BOUNDARY 2 /*Numericalflux constructors and destructor*/ @@ -33,13 +31,13 @@ /* Intermediary */ - int j; - int pos1,pos2,pos3,pos4; - int num_nodes; + int pos1,pos2,pos3,pos4; + int num_nodes; /*numericalflux constructor data: */ - int numericalflux_elem_ids[2]; - int numericalflux_vertex_ids[2]; - int numericalflux_node_ids[4]; - int numericalflux_type; + int numericalflux_elem_ids[2]; + int numericalflux_vertex_ids[2]; + int numericalflux_node_ids[4]; + int numericalflux_type; + int numericalflux_degree; /*Get edge*/ @@ -58,5 +56,5 @@ else{ /* internal edge: connected to 2 elements */ - num_nodes=4; + num_nodes=4; numericalflux_type=InternalEnum; numericalflux_elem_ids[0]=e1; @@ -64,4 +62,7 @@ } + /*FIXME: hardcode element degree for now*/ + this->flux_degree= P1DGEnum; + /*1: Get vertices ids*/ numericalflux_vertex_ids[0]=i1; @@ -69,11 +70,8 @@ /*2: Get node ids*/ - if (numericalflux_type==InternalEnum){ - - /*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*/ + if(numericalflux_type==InternalEnum){ + /*Get the column where these ids are located in the index*/ pos1=pos2=pos3=pos4=UNDEF; - for(j=0;j<3;j++){ + for(int j=0;j<3;j++){ if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; @@ -83,5 +81,5 @@ _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF); - /*3: We have the id of the elements and the position of the vertices in the index + /* 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*(e1-1)+pos1; @@ -91,8 +89,7 @@ } else{ - - /*2: Get the column where these ids are located in the index*/ + /*Get the column where these ids are located in the index*/ pos1=pos2=UNDEF; - for(j=0;j<3;j++){ + for(int j=0;j<3;j++){ if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; @@ -100,5 +97,5 @@ _assert_(pos1!=UNDEF && pos2!=UNDEF); - /*3: We have the id of the elements and the position of the vertices in the index + /* 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*(e1-1)+pos1; @@ -106,17 +103,18 @@ } - /*Ok, we have everything to build the object: */ - this->id=numericalflux_id; - this->flux_type = numericalflux_type; + /*Assign object fields: */ + this->id = numericalflux_id; + this->flux_type = numericalflux_type; + this->flux_degree = numericalflux_degree; /*Hooks: */ - this->hnodes =new Hook(numericalflux_node_ids,num_nodes); - this->hvertices =new Hook(&numericalflux_vertex_ids[0],2); - this->helement =new Hook(numericalflux_elem_ids,1); // take only the first element for now - - //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. - this->parameters=NULL; - this->element=NULL; - this->nodes=NULL; + this->hnodes = new Hook(numericalflux_node_ids,num_nodes); + this->hvertices = new Hook(&numericalflux_vertex_ids[0],2); + this->helement = new Hook(numericalflux_elem_ids,1); // take only the first element for now + + /*other fields*/ + this->parameters = NULL; + this->element = NULL; + this->nodes = NULL; } /*}}}*/ @@ -139,4 +137,5 @@ numericalflux->id=this->id; numericalflux->flux_type=this->flux_type; + numericalflux->flux_degree=this->flux_degree; /*point parameters: */ @@ -161,4 +160,5 @@ _printf_(" id: " << id << "\n"); _printf_(" flux_type: " << this->flux_type<< "\n"); + _printf_(" flux_degree: " << this->flux_degree<< "\n"); hnodes->DeepEcho(); hvertices->DeepEcho(); @@ -175,4 +175,5 @@ _printf_(" id: " << id << "\n"); _printf_(" flux_type: " << this->flux_type<< "\n"); + _printf_(" flux_degree: " << this->flux_degree<< "\n"); hnodes->Echo(); hvertices->Echo(); @@ -193,4 +194,5 @@ MARSHALLING(id); MARSHALLING(flux_type); + MARSHALLING(flux_degree); if(marshall_direction==MARSHALLING_BACKWARD){ @@ -300,14 +302,6 @@ _assert_(nodes); - switch(this->flux_type){ - case InternalEnum: - for(int i=0;iLid(); - return; - case BoundaryEnum: - for(int i=0;iLid(); - return; - default: - _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); - } + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iLid(); } /*}}}*/ @@ -317,25 +311,32 @@ _assert_(nodes); - switch(this->flux_type){ - case InternalEnum: - for(int i=0;iSid(); - return; - case BoundaryEnum: - for(int i=0;iSid(); - return; - default: - _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); - } + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iSid(); } /*}}}*/ int Numericalflux::GetNumberOfNodes(void){/*{{{*/ - switch(this->flux_type){ - case InternalEnum: - return NUMNODES_INTERNAL; - case BoundaryEnum: - return NUMNODES_BOUNDARY; - default: - _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); + if(this->flux_degree==P0DGEnum){ + switch(this->flux_type){ + case InternalEnum: + return 2; + case BoundaryEnum: + return 1; + default: + _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); + } + } + else if(this->flux_degree==P1DGEnum){ + switch(this->flux_type){ + case InternalEnum: + return 4; + case BoundaryEnum: + return 2; + default: + _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); + } + } + else{ + _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet"); } @@ -474,20 +475,12 @@ ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/ - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble L[numdof]; - IssmDouble Ke_g[numdof][numdof]; - GaussTria *gauss; - - /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke = NULL; - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; + IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; /*Retrieve all inputs and parameters*/ @@ -498,8 +491,8 @@ /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + + GaussTria* gauss=new GaussTria(); gauss->GaussEdgeCenter(index1,index2); vxaverage_input->GetInputValue(&mean_vx,gauss); @@ -507,19 +500,21 @@ delete gauss; - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN<=0){ + IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if(UdotN<=0){ return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ } - else{ - Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters); - } + + /*Initialize Element vector and other vectors*/ + int numnodes = this->GetNumberOfNodes(); + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *basis = xNew(numnodes); /* Start looping on the number of gaussian points: */ gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ + for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); vxaverage_input->GetInputValue(&vx,gauss); @@ -529,13 +524,13 @@ DL=gauss->weight*Jdet*UdotN; - TripleMultiply(&L[0],1,numdof,1, - &DL,1,1,0, - &L[0],1,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + for(int i=0;ivalues[i*numnodes+j]+=DL*basis[i]*basis[j]; + } + } } /*Clean up and return*/ + xDelete(basis); delete gauss; return Ke; @@ -543,23 +538,21 @@ /*}}}*/ ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/ - - /* constants*/ - const int numdof=NDOF1*NUMNODES_INTERNAL; - - /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL1,DL2,Jdet,vx,vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble B[numdof]; - IssmDouble Bprime[numdof]; - IssmDouble Ke_g1[numdof][numdof]; - IssmDouble Ke_g2[numdof][numdof]; - GaussTria *gauss; /*Initialize Element matrix and return if necessary*/ Tria* tria=(Tria*)element; if(!tria->IsIceInElement()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters); + + /* Intermediaries*/ + IssmDouble DL1,DL2,Jdet,vx,vy,UdotN; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + + /*Fetch number of nodes for this flux*/ + int numnodes = this->GetNumberOfNodes(); + + /*Initialize variables*/ + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *B = xNew(numnodes); + IssmDouble *Bprime = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -570,8 +563,8 @@ /* Start looping on the number of gaussian points: */ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + GaussTria* gauss=new GaussTria(index1,index2,2); + for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); @@ -587,18 +580,15 @@ DL2=gauss->weight*Jdet*fabs(UdotN)/2; - TripleMultiply(&B[0],1,numdof,1, - &DL1,1,1,0, - &Bprime[0],1,numdof,0, - &Ke_g1[0][0],0); - TripleMultiply(&B[0],1,numdof,1, - &DL2,1,1,0, - &B[0],1,numdof,0, - &Ke_g2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g1[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_g2[i][j]; + for(int i=0;ivalues[i*numnodes+j]+=DL1*B[i]*Bprime[j]; + Ke->values[i*numnodes+j]+=DL2*B[i]*B[j]; + } + } } /*Clean up and return*/ + xDelete(B); + xDelete(Bprime); delete gauss; return Ke; @@ -619,24 +609,16 @@ ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/ - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; - - /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble L[numdof]; - IssmDouble Ke_g[numdof][numdof]; - GaussTria *gauss; - /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke = NULL; Tria* tria=(Tria*)element; if(!tria->IsIceInElement()) return NULL; + /* Intermediaries*/ + IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&dt,TimesteppingTimeStepEnum); + IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); @@ -644,8 +626,8 @@ /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + + GaussTria* gauss=new GaussTria(); gauss->GaussEdgeCenter(index1,index2); vxaverage_input->GetInputValue(&mean_vx,gauss); @@ -653,19 +635,21 @@ delete gauss; - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN<=0){ + IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if(UdotN<=0){ return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ } - else{ - Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters); - } + + /*Initialize Element vector and other vectors*/ + int numnodes = this->GetNumberOfNodes(); + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *basis = xNew(numnodes); /* Start looping on the number of gaussian points: */ gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ + for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); vxaverage_input->GetInputValue(&vx,gauss); @@ -675,13 +659,13 @@ DL=gauss->weight*Jdet*dt*UdotN; - TripleMultiply(&L[0],1,numdof,1, - &DL,1,1,0, - &L[0],1,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + for(int i=0;ivalues[i*numnodes+j]+=DL*basis[i]*basis[j]; + } + } } /*Clean up and return*/ + xDelete(basis); delete gauss; return Ke; @@ -689,27 +673,25 @@ /*}}}*/ ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/ - - /* constants*/ - const int numdof=NDOF1*NUMNODES_INTERNAL; - - /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL1,DL2,Jdet,dt,vx,vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble B[numdof]; - IssmDouble Bprime[numdof]; - IssmDouble Ke_g1[numdof][numdof]; - IssmDouble Ke_g2[numdof][numdof]; - GaussTria *gauss; /*Initialize Element matrix and return if necessary*/ Tria* tria=(Tria*)element; if(!tria->IsIceInElement()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters); + + /* Intermediaries*/ + IssmDouble DL1,DL2,Jdet,vx,vy,UdotN; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + + /*Fetch number of nodes for this flux*/ + int numnodes = this->GetNumberOfNodes(); + + /*Initialize variables*/ + ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters); + IssmDouble *B = xNew(numnodes); + IssmDouble *Bprime = xNew(numnodes); /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&dt,TimesteppingTimeStepEnum); + IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); Input* vxaverage_input=tria->inputs->GetInput(VxEnum); Input* vyaverage_input=tria->inputs->GetInput(VyEnum); @@ -717,8 +699,8 @@ /* Start looping on the number of gaussian points: */ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + GaussTria* gauss=new GaussTria(index1,index2,2); + for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); @@ -734,18 +716,15 @@ DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2; - TripleMultiply(&B[0],1,numdof,1, - &DL1,1,1,0, - &Bprime[0],1,numdof,0, - &Ke_g1[0][0],0); - TripleMultiply(&B[0],1,numdof,1, - &DL2,1,1,0, - &B[0],1,numdof,0, - &Ke_g2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g1[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_g2[i][j]; + for(int i=0;ivalues[i*numnodes+j]+=DL1*B[i]*Bprime[j]; + Ke->values[i*numnodes+j]+=DL2*B[i]*B[j]; + } + } } /*Clean up and return*/ + xDelete(B); + xDelete(Bprime); delete gauss; return Ke; @@ -772,51 +751,45 @@ ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/ - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; + /*Initialize Load Vector and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; /* Intermediaries*/ - int i,ig,index1,index2; - IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness; + IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness; IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble normal[2]; - IssmDouble L[numdof]; - GaussTria *gauss; - - /*Initialize Load Vector and return if necessary*/ - ElementVector* pe = NULL; - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); - Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); - Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vxaverage_input = tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); + Input* vyaverage_input = tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); + Input* thickness_input = tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input); GetNormal(&normal[0],xyz_list); /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + GaussTria* gauss=new GaussTria(); gauss->GaussEdgeCenter(index1,index2); vxaverage_input->GetInputValue(&mean_vx,gauss); vyaverage_input->GetInputValue(&mean_vy,gauss); delete gauss; - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN>0){ + IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if(UdotN>0){ return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ } - else{ - pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters); - } + + /*Initialize Load Vector */ + int numnodes = this->GetNumberOfNodes(); + ElementVector *pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble *basis = xNew(numnodes); /* Start looping on the number of gaussian points: */ gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ + for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); vxaverage_input->GetInputValue(&vx,gauss); @@ -828,8 +801,9 @@ DL= - gauss->weight*Jdet*UdotN*thickness; - for(i=0;ivalues[i] += DL*L[i]; + for(int i=0;ivalues[i] += DL*basis[i]; } /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; @@ -857,53 +831,46 @@ ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/ - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; + /*Initialize Load Vector and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; /* Intermediaries*/ - int i,ig,index1,index2; - IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble L[numdof]; - GaussTria *gauss; - - /*Initialize Load Vector and return if necessary*/ - ElementVector* pe = NULL; - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; + IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vxaverage_input =tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); - Input* vyaverage_input =tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); - Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input); + IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum); + Input* vxaverage_input = tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); + Input* vyaverage_input = tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); + Input* spcthickness_input = tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input); GetNormal(&normal[0],xyz_list); /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); + GaussTria* gauss=new GaussTria(); gauss->GaussEdgeCenter(index1,index2); vxaverage_input->GetInputValue(&mean_vx,gauss); vyaverage_input->GetInputValue(&mean_vy,gauss); delete gauss; - - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN>0){ + IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if(UdotN>0){ return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ } - else{ - pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters); - } + + /*Initialize Load Vector */ + int numnodes = this->GetNumberOfNodes(); + ElementVector *pe = new ElementVector(nodes,numnodes,this->parameters); + IssmDouble *basis = xNew(numnodes); /* Start looping on the number of gaussian points: */ gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ + for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement()); vxaverage_input->GetInputValue(&vx,gauss); @@ -916,8 +883,9 @@ DL= - gauss->weight*Jdet*dt*UdotN*thickness; - for(i=0;ivalues[i] += DL*L[i]; + for(int i=0;ivalues[i] += DL*basis[i]; } /*Clean up and return*/ + xDelete(basis); delete gauss; return pe; @@ -935,10 +903,9 @@ /*Build unit outward pointing vector*/ IssmDouble vector[2]; - IssmDouble norm; vector[0]=xyz_list[1][0] - xyz_list[0][0]; vector[1]=xyz_list[1][1] - xyz_list[0][1]; - norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); + IssmDouble norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); normal[0]= + vector[1]/norm; Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 24087) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 24088) @@ -21,4 +21,5 @@ int id; int flux_type; + int flux_degree; /*Hooks*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 24087) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 24088) @@ -50,4 +50,7 @@ if(iomodel->meshelementtype==TriaEnum){ switch(finite_element){ + case P0DGEnum: + numnodes = iomodel->numberofelements; + break; case P1Enum: numnodes = iomodel->numberofvertices; @@ -173,4 +176,8 @@ if(iomodel->meshelementtype==TriaEnum){ switch(finite_element){ + case P0DGEnum: + element_numnodes=1; + element_node_ids[0]=i; + break; case P1Enum: element_numnodes=3; Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24087) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24088) @@ -1135,4 +1135,5 @@ syn keyword cConstant OptionEnum syn keyword cConstant P0ArrayEnum +syn keyword cConstant P0DGEnum syn keyword cConstant P1DGEnum syn keyword cConstant P1P1Enum @@ -1280,5 +1281,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1287,6 +1287,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1301,6 +1301,6 @@ syn keyword cType ElementHook syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1309,10 +1309,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1322,5 +1321,4 @@ syn keyword cType Input syn keyword cType Inputs -syn keyword cType IntArrayInput syn keyword cType IntInput syn keyword cType IntMatParam @@ -1330,5 +1328,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1341,5 +1338,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1354,6 +1350,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1368,10 +1364,10 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType Seg syn keyword cType SegInput +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1393,4 +1389,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1411,6 +1411,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GLheightadvectionAnalysis syn keyword cType GiaIvinsAnalysis -syn keyword cType GLheightadvectionAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24087) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24088) @@ -1133,4 +1133,5 @@ OptionEnum, P0ArrayEnum, + P0DGEnum, P1DGEnum, P1P1Enum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24087) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24088) @@ -1137,4 +1137,5 @@ case OptionEnum : return "Option"; case P0ArrayEnum : return "P0Array"; + case P0DGEnum : return "P0DG"; case P1DGEnum : return "P1DG"; case P1P1Enum : return "P1P1"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24087) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24088) @@ -1164,4 +1164,5 @@ else if (strcmp(name,"Option")==0) return OptionEnum; else if (strcmp(name,"P0Array")==0) return P0ArrayEnum; + else if (strcmp(name,"P0DG")==0) return P0DGEnum; else if (strcmp(name,"P1DG")==0) return P1DGEnum; else if (strcmp(name,"P1P1")==0) return P1P1Enum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"StringParam")==0) return StringParamEnum; else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum; - else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; + if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum; + else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum; else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum; else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;