Changeset 4422
- Timestamp:
- 07/07/10 08:42:38 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r4421 r4422 4714 4714 double Tria::GetAreaCoordinate(double x, double y, int which_one){ 4715 4715 4716 double area=0; 4717 const int numgrids=3; 4718 double xyz_list[numgrids][3]; 4719 double x1,y1,x2,y2,x3,y3; 4716 /*Intermediaries*/ 4717 double area = 0; 4718 const int numgrids = 3; 4719 double xyz_list[numgrids][3]; 4720 double x1,y1,x2,y2,x3,y3; 4720 4721 4721 4722 /*Get area: */ … … 4728 4729 x3=xyz_list[2][0]; y3=xyz_list[2][1]; 4729 4730 4730 if(which_one==1){4731 /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/4732 return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area;4733 }4734 else if(which_one==2){4735 /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/4736 return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;4737 }4738 else if(which_one==3){4739 /*Get third area coordinate 1-area1-area2: */4740 return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;4741 }4742 else ISSMERROR("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!");4731 switch(which_one){ 4732 case 1: 4733 /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ 4734 return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; 4735 case 2: 4736 /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ 4737 return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; 4738 case 3: 4739 /*Get third area coordinate 1-area1-area2: */ 4740 return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; 4741 default: 4742 ISSMERROR("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"); 4743 } 4743 4744 } 4744 4745 /*}}}*/ … … 4953 4954 y3=*(xyz_list+3*2+1); 4954 4955 4955 4956 4956 *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); 4957 4958 4957 4959 4958 if(Jdet<0){ … … 5149 5148 } 5150 5149 /*}}}*/ 5151 /*FUNCTION Tria::GetParameterValue {{{1*/5150 /*FUNCTION Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3) {{{1*/ 5152 5151 void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){ 5153 5152 … … 5167 5166 /*Assign output pointers:*/ 5168 5167 *pp=p; 5168 } 5169 /*}}}*/ 5170 /*FUNCTION Tria::GetParameterValue( {{{1*/ 5171 void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype){ 5172 5173 /*Output*/ 5174 double value; 5175 5176 /*Intermediaries*/ 5177 const int numnodes=3; 5178 int grid1=-1,grid2=-1; 5179 int grid3; 5180 int i; 5181 double gauss_tria[numnodes]; 5182 5183 /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */ 5184 for(i=0;i<numnodes;i++){ 5185 if (node1==nodes[i]) grid1=i; 5186 if (node2==nodes[i]) grid2=i; 5187 } 5188 5189 /*Reverse grid1 and 2 if necessary*/ 5190 if (grid1>grid2){ 5191 5192 /*Reverse grid1 and grid2*/ 5193 grid3=grid1; grid1=grid2; grid2=grid3; 5194 5195 /*Change segment gauss point (between -1 and +1)*/ 5196 gauss_seg=-gauss_seg; 5197 } 5198 5199 /*Build Triangle Gauss point*/ 5200 if (grid1==0 && grid2==1){ 5201 /*gauss_seg is between 0 and 1*/ 5202 gauss_tria[0]=0.5*(1.-gauss_seg); 5203 gauss_tria[1]=1.-0.5*(1.-gauss_seg); 5204 gauss_tria[2]=0.; 5205 } 5206 else if (grid1==0 && grid2==2){ 5207 /*gauss_seg is between 0 and 2*/ 5208 gauss_tria[0]=0.5*(1.-gauss_seg); 5209 gauss_tria[1]=0.; 5210 gauss_tria[2]=1.-0.5*(1.-gauss_seg); 5211 } 5212 else if (grid1==1 && grid2==2){ 5213 /*gauss_seg is between 1 and 2*/ 5214 gauss_tria[0]=0.; 5215 gauss_tria[1]=0.5*(1.-gauss_seg); 5216 gauss_tria[2]=1.-0.5*(1.-gauss_seg); 5217 } 5218 else{ 5219 ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes"); 5220 } 5221 5222 /*OK, now we can call input method*/ 5223 this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype); 5224 5225 /*Assign output pointers:*/ 5226 *pvalue=value; 5169 5227 } 5170 5228 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Tria.h
r4402 r4422 150 150 void GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3); 151 151 void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); 152 void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype); 152 153 void GetSolutionFromInputsDiagnosticHoriz(Vec solution); 153 154 void GetSolutionFromInputsAdjoint(Vec solution); -
issm/trunk/src/c/objects/FemModel.cpp
r4402 r4422 143 143 144 144 /*Numerics: */ 145 /*FUNCTION FemModel::SetCurrentConfiguration {{{1*/145 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{1*/ 146 146 void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){ 147 147 … … 174 174 } 175 175 /*}}}1*/ 176 /*FUNCTION FemModel::SetCurrentConfiguration {{{1*/176 /*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{1*/ 177 177 void FemModel::SetCurrentConfiguration(int configuration_type){ 178 178 -
issm/trunk/src/c/objects/Loads/Numericalflux.cpp
r4405 r4422 432 432 433 433 //Get vx, vy and v.n 434 t rias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum);435 t rias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum);434 this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); 435 this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); 436 436 437 437 UdotN=vx*normal[0]+vy*normal[1]; … … 536 536 GetNormal(&normal[0],xyz_list); 537 537 538 /*Check wether it is an inflow or outflow BC */539 t rias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum);540 t rias[0]->inputs->GetParameterValue(&mean_vy, nodes[0],nodes[1],.5,VyAverageEnum);538 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ 539 this->GetParameterValue(&mean_vx,0.,VxAverageEnum); 540 this->GetParameterValue(&mean_vy,0.,VyAverageEnum); 541 541 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 542 542 if (UdotN<=0){ … … 560 560 561 561 //Get vx, vy and v.n 562 t rias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum);563 t rias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum);562 this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); 563 this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); 564 564 565 565 UdotN=vx*normal[0]+vy*normal[1]; … … 665 665 GetNormal(&normal[0],xyz_list); 666 666 667 /*Check wether it is an inflow or outflow BC */668 t rias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum);669 t rias[0]->inputs->GetParameterValue(&mean_vy, nodes[0],nodes[1],.5,VyAverageEnum);667 /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ 668 this->GetParameterValue(&mean_vx,0.,VxAverageEnum); 669 this->GetParameterValue(&mean_vy,0.,VyAverageEnum); 670 670 UdotN=mean_vx*normal[0]+mean_vy*normal[1]; 671 671 if (UdotN>0){ … … 689 689 690 690 //Get vx, vy and v.n 691 t rias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum);692 t rias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum);693 t rias[0]->inputs->GetParameterValue(&thickness, nodes[0],nodes[1],gauss_coord,ThicknessEnum);691 this->GetParameterValue(&vx,gauss_coord,VxAverageEnum); 692 this->GetParameterValue(&vy,gauss_coord,VyAverageEnum); 693 this->GetParameterValue(&thickness,gauss_coord,ThicknessEnum); 694 694 695 695 UdotN=vx*normal[0]+vy*normal[1]; … … 836 836 /*}}}*/ 837 837 /*FUNCTION Numericalflux::GetParameterValue {{{1*/ 838 void Numericalflux::GetParameterValue(double* pp, double* plist, double gauss_coord){ 839 840 /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian 841 * point specifie by gauss_l1l2l3: */ 842 843 /*nodal functions: */ 844 double l1l4[4]; 838 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype){ 845 839 846 840 /*output: */ 847 double p; 848 849 GetNodalFunctions(&l1l4[0],gauss_coord); 850 851 p=l1l4[0]*plist[0]+l1l4[1]*plist[1]; 841 double value1,value2,value; 842 int type; 843 844 /*dynamic objects pointed to by hooks: */ 845 Tria** trias=NULL; 846 Node** nodes=NULL; 847 848 /*recover type: */ 849 inputs->GetParameterValue(&type,TypeEnum); 850 851 852 /*recover objects from hooks: */ 853 trias=(Tria**)helements->deliverp(); 854 nodes=(Node**)hnodes->deliverp(); 855 856 if (type==InternalEnum){ 857 /*Get value on Element 1 and 2 and take the average*/ 858 trias[0]->GetParameterValue(&value1,nodes[0],nodes[1],gauss_coord,enumtype); 859 trias[1]->GetParameterValue(&value2,nodes[2],nodes[3],gauss_coord,enumtype); 860 value=0.5*(value1+value2); 861 } 862 else{ 863 /*Get value on Element*/ 864 trias[0]->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,enumtype); 865 } 852 866 853 867 /*Assign output pointers:*/ 854 *p p=p;855 } 856 /*}}}*/ 868 *pvalue=value; 869 } 870 /*}}}*/ -
issm/trunk/src/c/objects/Loads/Numericalflux.h
r4396 r4422 68 68 void GetDofList(int* doflist,int* pnumberofdofs); 69 69 void GetNormal(double* normal,double xyz_list[4][3]); 70 void GetParameterValue(double* p p, double* plist, double gauss_coord);70 void GetParameterValue(double* pvalue, double gauss_coord,int enumtype); 71 71 void CreateKMatrixInternal(Mat Kgg); 72 72 void CreateKMatrixBoundary(Mat Kgg);
Note:
See TracChangeset
for help on using the changeset viewer.