Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4903) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4904) @@ -1962,4 +1962,8 @@ if(!onbed)return; + /*Depth Averaging Vx and Vy*/ + this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); + this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); + /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. @@ -1990,4 +1994,8 @@ /*Is this element on the bed? :*/ if(!onbed)return; + + /*Depth Averaging Vx and Vy*/ + this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); + this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); /*Spawn Tria element from the base of the Penta: */ @@ -2691,5 +2699,4 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrixPrognostic {{{1*/ - void Penta::CreateKMatrixPrognostic(Mat Kgg){ @@ -2710,4 +2717,8 @@ /*Is this element on the bed? :*/ if(!onbed)return; + + /*Depth Averaging Vx and Vy*/ + this->InputDepthAverageAtBase(VxEnum,VxAverageEnum); + this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); /*Spawn Tria element from the base of the Penta: */ @@ -2720,5 +2731,4 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrixSlope {{{1*/ - void Penta::CreateKMatrixSlope(Mat Kgg){ @@ -3692,5 +3702,4 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorPrognostic {{{1*/ - void Penta::CreatePVectorPrognostic( Vec pg){ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4903) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4904) @@ -2321,9 +2321,9 @@ double v_gauss[2]={0.0}; - double K[2][2]={0.0}; double KDL[2][2]={0.0}; int dofs[2]={0,1}; int found=0; + int dim; /*inputs: */ @@ -2338,8 +2338,15 @@ /*retrieve some parameters: */ this->parameters->FindParam(&artdiff,ArtDiffEnum); + this->parameters->FindParam(&dim,DimEnum); /*Retrieve all inputs we will be needed*/ - vxaverage_input=inputs->GetInput(VxAverageEnum); - vyaverage_input=inputs->GetInput(VyAverageEnum); + if(dim==2){ + vxaverage_input=inputs->GetInput(VxEnum); + vyaverage_input=inputs->GetInput(VyEnum); + } + else{ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); + } //Create Artificial diffusivity once for all if requested @@ -2476,4 +2483,5 @@ int dofs[1]={0}; int found; + int dim; /*inputs: */ @@ -2484,4 +2492,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); + this->parameters->FindParam(&dim,DimEnum); /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ @@ -2489,6 +2498,8 @@ /*Retrieve all inputs we will be needed*/ - vx_input=inputs->GetInput(VxEnum); - vy_input=inputs->GetInput(VyEnum); + if(dim==2){ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + } /* Start looping on the number of gaussian points: */ @@ -2589,12 +2600,15 @@ int dofs[2]={0,1}; int found=0; + int dim; /*inputs: */ Input* vxaverage_input=NULL; Input* vyaverage_input=NULL; + Input* surface_input=NULL; bool artdiff; /*retrieve some parameters: */ this->parameters->FindParam(&artdiff,ArtDiffEnum); + this->parameters->FindParam(&dim,DimEnum); /* Get node coordinates and dof list: */ @@ -2602,11 +2616,22 @@ GetDofList(&doflist[0],&numberofdofspernode); + /*Retrieve all inputs we will be needed*/ + surface_input=inputs->GetInput(SurfaceEnum); + if(dim==2){ + vxaverage_input=inputs->GetInput(VxEnum); + vyaverage_input=inputs->GetInput(VyEnum); + } + else{ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); + } + /*Modify z so that it reflects the surface*/ - inputs->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3,SurfaceEnum); + surface_input->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3); for(i=0;iGetParameterAverage(&nx,VxAverageEnum); - inputs->GetParameterAverage(&ny,VyAverageEnum); + vxaverage_input->GetParameterAverage(&nx); + vyaverage_input->GetParameterAverage(&ny); if(nx==0 && ny==0){ SurfaceNormal(&surface_normal[0],xyz_list); @@ -2621,8 +2646,4 @@ nx=nx/norm; ny=ny/norm; - - /*Retrieve all inputs we will be needed*/ - vxaverage_input=inputs->GetInput(VxAverageEnum); - vyaverage_input=inputs->GetInput(VyAverageEnum); //Create Artificial diffusivity once for all if requested @@ -2704,11 +2725,8 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - -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); - } /*}}}*/ @@ -3284,4 +3302,5 @@ int dofs[2]={0,1}; int found; + int dim; /*inputs: */ @@ -3293,4 +3312,5 @@ /*retrieve some parameters: */ this->parameters->FindParam(&dt,DtEnum); + this->parameters->FindParam(&dim,DimEnum); this->parameters->FindParam(&artdiff,ArtDiffEnum); @@ -3300,6 +3320,12 @@ /*Retrieve all inputs we will be needing: */ - vxaverage_input=inputs->GetInput(VxAverageEnum); - vyaverage_input=inputs->GetInput(VyAverageEnum); + if(dim==2){ + vxaverage_input=inputs->GetInput(VxEnum); + vyaverage_input=inputs->GetInput(VyEnum); + } + else{ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); + } //Create Artificial diffusivity once for all if requested @@ -3405,6 +3431,4 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - -cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -3453,4 +3477,5 @@ int dofs[1]={0}; int found; + int dim; /*inputs: */ @@ -3461,4 +3486,5 @@ /*retrieve some parameters: */ this->parameters->FindParam(&dt,DtEnum); + this->parameters->FindParam(&dim,DimEnum); /* Get node coordinates and dof list: */ @@ -3470,6 +3496,12 @@ /*Retrieve all inputs we will be needed*/ - vxaverage_input=inputs->GetInput(VxAverageEnum); - vyaverage_input=inputs->GetInput(VyAverageEnum); + if(dim==2){ + vxaverage_input=inputs->GetInput(VxEnum); + vyaverage_input=inputs->GetInput(VyEnum); + } + else{ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); + } /* Start looping on the number of gaussian points: */ @@ -3525,6 +3557,4 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - -cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -3710,5 +3740,4 @@ void Tria::CreatePVectorBalancedthickness_CG(Vec pg ){ - /* local declarations */ int i,j; @@ -3781,10 +3810,8 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); -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); - } /*}}}*/ @@ -3866,15 +3893,12 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); -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::CreatePVectorBalancedvelocities {{{1*/ void Tria::CreatePVectorBalancedvelocities(Vec pg){ - /* local declarations */ @@ -4814,5 +4838,4 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); -cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -4901,5 +4924,4 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); -cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -5470,4 +5492,64 @@ /*OK, now we can call input method*/ this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype); + + /*Assign output pointers:*/ + *pvalue=value; +} +/*}}}*/ +/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/ +void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){ + + /*Output*/ + double value; + + /*Intermediaries*/ + const int numnodes=3; + int grid1=-1,grid2=-1; + int grid3; + int i; + double gauss_tria[numnodes]; + + /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */ + ISSMASSERT(nodes); + for(i=0;igrid2){ + + /*Reverse grid1 and grid2*/ + grid3=grid1; grid1=grid2; grid2=grid3; + + /*Change segment gauss point (between -1 and +1)*/ + gauss_seg=-gauss_seg; + } + + /*Build Triangle Gauss point*/ + if (grid1==0 && grid2==1){ + /*gauss_seg is between 0 and 1*/ + gauss_tria[0]=0.5*(1.-gauss_seg); + gauss_tria[1]=1.-0.5*(1.-gauss_seg); + gauss_tria[2]=0.; + } + else if (grid1==0 && grid2==2){ + /*gauss_seg is between 0 and 2*/ + gauss_tria[0]=0.5*(1.-gauss_seg); + gauss_tria[1]=0.; + gauss_tria[2]=1.-0.5*(1.-gauss_seg); + } + else if (grid1==1 && grid2==2){ + /*gauss_seg is between 1 and 2*/ + gauss_tria[0]=0.; + gauss_tria[1]=0.5*(1.-gauss_seg); + gauss_tria[2]=1.-0.5*(1.-gauss_seg); + } + else{ + ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes"); + } + + /*OK, now we can call input method*/ + input_in->GetParameterValue(&value, &gauss_tria[0]); /*Assign output pointers:*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4903) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4904) @@ -145,4 +145,5 @@ void GetParameterValue(double* pvalue,Node* node,int enumtype); void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype); + void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in); void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsAdjointHoriz(Vec solution); Index: /issm/trunk/src/c/objects/Inputs/Input.h =================================================================== --- /issm/trunk/src/c/objects/Inputs/Input.h (revision 4903) +++ /issm/trunk/src/c/objects/Inputs/Input.h (revision 4904) @@ -12,5 +12,4 @@ class Node; class ElementResult; -#include "../Node.h" /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 4903) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 4904) @@ -401,13 +401,19 @@ double dt; int found; + int dim; /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; + Tria* tria=NULL; + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; /*Retrieve parameters: */ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&dim,DimEnum); /*recover objects from hooks: */ nodes=(Node**)hnodes->deliverp(); + tria=(Tria*)helement->delivers(); /*recover parameters: */ @@ -428,4 +434,10 @@ 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); + } + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ num_gauss=2; @@ -443,6 +455,6 @@ //Get vx, vy and v.n - this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); - this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); + this->GetParameterValue(&vx,gauss_coord,vxaverage_input); + this->GetParameterValue(&vy,gauss_coord,vyaverage_input); UdotN=vx*normal[0]+vy*normal[1]; @@ -518,13 +530,19 @@ int dofs[1]={0}; int found; + int dim; /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; + Tria* tria=NULL; + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; /*Retrieve parameters: */ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&dim,DimEnum); /*recover objects from hooks: */ nodes=(Node**)hnodes->deliverp(); + tria=(Tria*)helement->delivers(); /*recover parameters: */ @@ -545,7 +563,13 @@ 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); + } + /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - this->GetParameterValue(&mean_vx,0.,VxAverageEnum); - this->GetParameterValue(&mean_vy,0.,VyAverageEnum); + this->GetParameterValue(&mean_vx,0.,vxaverage_input); + this->GetParameterValue(&mean_vy,0.,vyaverage_input); UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN<=0){ @@ -569,6 +593,6 @@ //Get vx, vy and v.n - this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); - this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); + this->GetParameterValue(&vx,gauss_coord,vxaverage_input); + this->GetParameterValue(&vy,gauss_coord,vyaverage_input); UdotN=vx*normal[0]+vy*normal[1]; @@ -645,13 +669,20 @@ int dofs[1]={0}; int found; + int dim; /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; + Tria* tria=NULL; + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; + Input* thickness_input=NULL; /*recover objects from hooks: */ nodes=(Node**)hnodes->deliverp(); + tria=(Tria*)helement->delivers(); /*Retrieve parameters: */ this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&dim,DimEnum); /*recover parameters: */ @@ -667,4 +698,11 @@ } + /*Retrieve all inputs we will be needing: */ + if(dim==2){ + vxaverage_input=tria->inputs->GetInput(VxEnum); + vyaverage_input=tria->inputs->GetInput(VyEnum); + } + thickness_input=tria->inputs->GetInput(ThicknessEnum); + /* Get node coordinates, dof list and normal vector: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -673,6 +711,6 @@ /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - this->GetParameterValue(&mean_vx,0.,VxAverageEnum); - this->GetParameterValue(&mean_vy,0.,VyAverageEnum); + this->GetParameterValue(&mean_vx,0.,vxaverage_input); + this->GetParameterValue(&mean_vy,0.,vyaverage_input); UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN>0){ @@ -696,7 +734,7 @@ //Get vx, vy and v.n - this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); - this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); - this->GetParameterValue(&thickness,gauss_coord,ThicknessEnum); + this->GetParameterValue(&vx,gauss_coord,vxaverage_input); + this->GetParameterValue(&vy,gauss_coord,vyaverage_input); + this->GetParameterValue(&thickness,gauss_coord,thickness_input); UdotN=vx*normal[0]+vy*normal[1]; @@ -846,6 +884,5 @@ /*output: */ - double value1,value2,value; - int type; + double value; /*dynamic objects pointed to by hooks: */ @@ -853,8 +890,4 @@ Node** nodes=NULL; - /*recover type: */ - inputs->GetParameterValue(&type,TypeEnum); - - /*recover objects from hooks: */ tria=(Tria*)helement->delivers(); @@ -868,2 +901,23 @@ } /*}}}*/ +/*FUNCTION Numericalflux::GetParameterValue {{{1*/ +void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){ + + /*output: */ + double value; + + /*dynamic objects pointed to by hooks: */ + Tria* tria=NULL; + Node** nodes=NULL; + + /*recover objects from hooks: */ + tria=(Tria*)helement->delivers(); + nodes=(Node**)hnodes->deliverp(); + + /*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 4903) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 4904) @@ -70,4 +70,5 @@ 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);