Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4421) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4422) @@ -4714,8 +4714,9 @@ double Tria::GetAreaCoordinate(double x, double y, int which_one){ - double area=0; - const int numgrids=3; - double xyz_list[numgrids][3]; - double x1,y1,x2,y2,x3,y3; + /*Intermediaries*/ + double area = 0; + const int numgrids = 3; + double xyz_list[numgrids][3]; + double x1,y1,x2,y2,x3,y3; /*Get area: */ @@ -4728,17 +4729,17 @@ x3=xyz_list[2][0]; y3=xyz_list[2][1]; - if(which_one==1){ - /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ - return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; - } - else if(which_one==2){ - /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ - return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; - } - else if(which_one==3){ - /*Get third area coordinate 1-area1-area2: */ - return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; - } - else ISSMERROR("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"); + switch(which_one){ + case 1: + /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ + return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; + case 2: + /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ + return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; + case 3: + /*Get third area coordinate 1-area1-area2: */ + return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; + default: + ISSMERROR("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"); + } } /*}}}*/ @@ -4953,7 +4954,5 @@ y3=*(xyz_list+3*2+1); - *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); - if(Jdet<0){ @@ -5149,5 +5148,5 @@ } /*}}}*/ -/*FUNCTION Tria::GetParameterValue {{{1*/ +/*FUNCTION Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3) {{{1*/ void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){ @@ -5167,4 +5166,63 @@ /*Assign output pointers:*/ *pp=p; +} +/*}}}*/ +/*FUNCTION Tria::GetParameterValue( {{{1*/ +void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype){ + + /*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: */ + 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*/ + this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype); + + /*Assign output pointers:*/ + *pvalue=value; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4421) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4422) @@ -150,4 +150,5 @@ void GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3); void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); + void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype); void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsAdjoint(Vec solution); Index: /issm/trunk/src/c/objects/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 4421) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 4422) @@ -143,5 +143,5 @@ /*Numerics: */ -/*FUNCTION FemModel::SetCurrentConfiguration{{{1*/ +/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{1*/ void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){ @@ -174,5 +174,5 @@ } /*}}}1*/ -/*FUNCTION FemModel::SetCurrentConfiguration{{{1*/ +/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{1*/ void FemModel::SetCurrentConfiguration(int configuration_type){ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 4421) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 4422) @@ -432,6 +432,6 @@ //Get vx, vy and v.n - trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum); - trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum); + this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); + this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); UdotN=vx*normal[0]+vy*normal[1]; @@ -536,7 +536,7 @@ GetNormal(&normal[0],xyz_list); - /*Check wether it is an inflow or outflow BC*/ - trias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum); - trias[0]->inputs->GetParameterValue(&mean_vy, nodes[0],nodes[1],.5,VyAverageEnum); + /*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); UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN<=0){ @@ -560,6 +560,6 @@ //Get vx, vy and v.n - trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum); - trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum); + this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); + this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); UdotN=vx*normal[0]+vy*normal[1]; @@ -665,7 +665,7 @@ GetNormal(&normal[0],xyz_list); - /*Check wether it is an inflow or outflow BC*/ - trias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum); - trias[0]->inputs->GetParameterValue(&mean_vy, nodes[0],nodes[1],.5,VyAverageEnum); + /*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); UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN>0){ @@ -689,7 +689,7 @@ //Get vx, vy and v.n - trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum); - trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum); - trias[0]->inputs->GetParameterValue(&thickness, nodes[0],nodes[1],gauss_coord,ThicknessEnum); + this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); + this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); + this->GetParameterValue(&thickness,gauss_coord,ThicknessEnum); UdotN=vx*normal[0]+vy*normal[1]; @@ -836,21 +836,35 @@ /*}}}*/ /*FUNCTION Numericalflux::GetParameterValue {{{1*/ -void Numericalflux::GetParameterValue(double* pp, double* plist, double gauss_coord){ - - /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian - * point specifie by gauss_l1l2l3: */ - - /*nodal functions: */ - double l1l4[4]; +void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype){ /*output: */ - double p; - - GetNodalFunctions(&l1l4[0],gauss_coord); - - p=l1l4[0]*plist[0]+l1l4[1]*plist[1]; + double value1,value2,value; + int type; + + /*dynamic objects pointed to by hooks: */ + Tria** trias=NULL; + Node** nodes=NULL; + + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + + /*recover objects from hooks: */ + trias=(Tria**)helements->deliverp(); + nodes=(Node**)hnodes->deliverp(); + + if (type==InternalEnum){ + /*Get value on Element 1 and 2 and take the average*/ + trias[0]->GetParameterValue(&value1,nodes[0],nodes[1],gauss_coord,enumtype); + trias[1]->GetParameterValue(&value2,nodes[2],nodes[3],gauss_coord,enumtype); + value=0.5*(value1+value2); + } + else{ + /*Get value on Element*/ + trias[0]->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,enumtype); + } /*Assign output pointers:*/ - *pp=p; -} -/*}}}*/ + *pvalue=value; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 4421) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 4422) @@ -68,5 +68,5 @@ void GetDofList(int* doflist,int* pnumberofdofs); void GetNormal(double* normal,double xyz_list[4][3]); - void GetParameterValue(double* pp, double* plist, double gauss_coord); + void GetParameterValue(double* pvalue, double gauss_coord,int enumtype); void CreateKMatrixInternal(Mat Kgg); void CreateKMatrixBoundary(Mat Kgg);