Index: /issm/trunk/src/c/objects/Gauss/GaussTria.cpp =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussTria.cpp (revision 5738) +++ /issm/trunk/src/c/objects/Gauss/GaussTria.cpp (revision 5739) @@ -56,5 +56,5 @@ weights=(double*)xmalloc(numgauss*sizeof(double)); - /*Reverse grid1 and 2 if necessary*/ + /*Reverse index1 and 2 if necessary*/ if (index1>index2){ index3=index1; index1=index2; index2=index3; @@ -150,4 +150,35 @@ coord2=ONETHIRD; coord3=ONETHIRD; + +} +/*}}}*/ +/*FUNCTION GaussTria::GaussEdgeCenter{{{1*/ +void GaussTria::GaussEdgeCenter(int index1,int index2){ + + int index3; + + /*Reverse index1 and 2 if necessary*/ + if (index1>index2){ + index3=index1; index1=index2; index2=index3; + } + + /*update static arrays*/ + if (index1==0 && index2==1){ + coord1=0.5; + coord2=0.5; + coord3=0.0; + } + else if (index1==0 && index2==2){ + coord1=0.5; + coord2=0.0; + coord3=0.5; + } + else if (index1==1 && index2==2){ + coord1=0.0; + coord2=0.5; + coord3=0.5; + } + else + ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2); } Index: /issm/trunk/src/c/objects/Gauss/GaussTria.h =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussTria.h (revision 5738) +++ /issm/trunk/src/c/objects/Gauss/GaussTria.h (revision 5739) @@ -42,4 +42,5 @@ void GaussVertex(int iv); void GaussCenter(void); + void GaussEdgeCenter(int index1,int index2); }; #endif /* _GAUSSTRIA_H_ */ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5738) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5739) @@ -20,5 +20,6 @@ /*Load macros*/ -#define NUMVERTICES 4 +#define NUMVERTICES_INTERNAL 4 +#define NUMVERTICES_BOUNDARY 2 #define NDOF1 1 @@ -400,20 +401,20 @@ /* constants*/ - const int numdof=NDOF1*NUMVERTICES; + const int numdof=NDOF1*NUMVERTICES_INTERNAL; /* Intermediaries*/ int i,j,ig,index1,index2,analysis_type; double DL1,DL2,Jdet,dt,vx,vy,UdotN; - double xyz_list[NUMVERTICES][3]; + double xyz_list[NUMVERTICES_INTERNAL][3]; double normal[2]; double B[numdof]; double Bprime[numdof]; - double Ke_gg1[numdof][numdof]; - double Ke_gg2[numdof][numdof]; - double Ke_gg[numdof][numdof] = {0.0}; + double Ke_g1[numdof][numdof]; + double Ke_g2[numdof][numdof]; + double Ke[numdof][numdof] = {0.0}; int *doflist = NULL; GaussTria *gauss; - GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES); + GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL); GetDofList(&doflist); @@ -454,15 +455,15 @@ &DL1,1,1,0, &Bprime[0],1,numdof,0, - &Ke_gg1[0][0],0); + &Ke_g1[0][0],0); TripleMultiply(&B[0],1,numdof,1, &DL2,1,1,0, &B[0],1,numdof,0, - &Ke_gg2[0][0],0); - - for(i=0;iinputs->GetInput(VxEnum); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&dim,DimEnum); - - /*recover objects from hooks: */ - Tria* tria=(Tria*)element; - - /*recover parameters: */ - if (analysis_type==PrognosticAnalysisEnum){ - parameters->FindParam(&dt,DtEnum); - } - else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){ - /*No transient term is involved*/ - dt=1; - } - else{ - ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); - } - - /* Get node coordinates and normal vector: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetNormal(&normal[0],xyz_list); - - /*Retrieve all inputs we will be needing: */ - if(dim==2){ - vxaverage_input=tria->inputs->GetInput(VxEnum); - vyaverage_input=tria->inputs->GetInput(VyEnum); + switch(analysis_type){ + case PrognosticAnalysisEnum: + parameters->FindParam(&dt,DtEnum); break; + case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + dt=1; break;/*No transient term is involved*/ + default: + ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); } /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - this->GetParameterValue(&mean_vx,0.,vxaverage_input); - this->GetParameterValue(&mean_vy,0.,vyaverage_input); + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + + gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + vxaverage_input->GetParameterValue(&mean_vx,gauss); + vyaverage_input->GetParameterValue(&mean_vy,gauss); + delete gauss; + UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN<=0){ @@ -551,53 +523,33 @@ } - /*Get dof list*/ GetDofList(&doflist); - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - num_gauss=2; - gaussSegment(&gauss_coords, &gauss_weights,num_gauss); - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetParameterValue(&vx,gauss_coord,vxaverage_input); - this->GetParameterValue(&vy,gauss_coord,vyaverage_input); - + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); + + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); UdotN=vx*normal[0]+vy*normal[1]; - - /*Get L*/ - GetL(&L[0],gauss_coord,numdof); - - /*Compute DLs*/ - DL=gauss_weight*Jdet*dt*UdotN; - - /*Do triple product*/ + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL=gauss->weight*Jdet*dt*UdotN; + TripleMultiply(&L[0],1,numdof,1, &DL,1,1,0, &L[0],1,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; iinputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input); + Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); + + /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/ + if (!thickness_input){ + thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + } + this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&dim,DimEnum); - - /*recover parameters: */ - if (analysis_type==PrognosticAnalysisEnum){ - parameters->FindParam(&dt,DtEnum); - } - else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){ - /*No transient term is involved*/ - dt=1; - } - else{ - ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); - } - - /*Retrieve all inputs we will be needing: */ - if(dim==2){ - vxaverage_input=tria->inputs->GetInput(VxEnum); - vyaverage_input=tria->inputs->GetInput(VyEnum); - } - /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/ - thickness_input=tria->inputs->GetInput(ThicknessObsEnum); - if (!thickness_input) thickness_input=tria->inputs->GetInput(ThicknessEnum); - - /* Get node coordinates and normal vector: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetNormal(&normal[0],xyz_list); + switch(analysis_type){ + case PrognosticAnalysisEnum: + parameters->FindParam(&dt,DtEnum); break; + case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + dt=1; break;/*No transient term is involved*/ + default: + ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type)); + } /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - this->GetParameterValue(&mean_vx,0.,vxaverage_input); - this->GetParameterValue(&mean_vy,0.,vyaverage_input); + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + + gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + vxaverage_input->GetParameterValue(&mean_vx,gauss); + vyaverage_input->GetParameterValue(&mean_vy,gauss); + delete gauss; UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN>0){ @@ -693,64 +617,30 @@ } - /*Get dof list*/ GetDofList(&doflist); - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - num_gauss=2; - gaussSegment(&gauss_coords, &gauss_weights,num_gauss); - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetParameterValue(&vx,gauss_coord,vxaverage_input); - this->GetParameterValue(&vy,gauss_coord,vyaverage_input); - this->GetParameterValue(&thickness,gauss_coord,thickness_input); - + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2); + + vxaverage_input->GetParameterValue(&vx,gauss); + vyaverage_input->GetParameterValue(&vy,gauss); + thickness_input->GetParameterValue(&thickness,gauss); UdotN=vx*normal[0]+vy*normal[1]; - - /*Get L*/ - GetL(&L[0],gauss_coord,numdof); - - /*Compute DL*/ - DL= - gauss_weight*Jdet*dt*UdotN*thickness; - - /* Add value into pe_g: */ - for( i=0; iGetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL= - gauss->weight*Jdet*dt*UdotN*thickness; + + for(i=0;iGetParameterValue(&value,nodes[0],nodes[1],gauss_coord,enumtype); - - /*Assign output pointers:*/ - *pvalue=value; -} -/*}}}*/ -/*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/ -void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){ - - /*output: */ - double value; - - /*recover objects from hooks: */ - Tria* tria=(Tria*)element; - - /*Get value on Element 1*/ - tria->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in); - - /*Assign output pointers:*/ - *pvalue=value; -} -/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 5738) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 5739) @@ -72,12 +72,6 @@ /*}}}*/ /*Numericalflux management:{{{1*/ - void GetJacobianDeterminant(double* pJdet,double xyz_list[4][3], double gauss_coord); - void GetNodalFunctions(double* l1l4, double gauss_coord); - void GetB(double* B, double gauss_coord); - void GetL(double* L, double gauss_coord,int numdof); void GetDofList(int** pdoflist); void GetNormal(double* normal,double xyz_list[4][3]); - void GetParameterValue(double* pvalue, double gauss_coord,int enumtype); - void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in); void CreateKMatrixInternal(Mat Kgg); void CreateKMatrixBoundary(Mat Kgg);