Changeset 4422


Ignore:
Timestamp:
07/07/10 08:42:38 (15 years ago)
Author:
Mathieu Morlighem
Message:

Better GetParameterValue for NumericalFlux

Location:
issm/trunk/src/c/objects
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r4421 r4422  
    47144714double Tria::GetAreaCoordinate(double x, double y, int which_one){
    47154715
    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;
    47204721
    47214722        /*Get area: */
    … …  
    47284729        x3=xyz_list[2][0]; y3=xyz_list[2][1];
    47294730
    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        }
    47434744}
    47444745/*}}}*/
    … …  
    49534954        y3=*(xyz_list+3*2+1);
    49544955
    4955 
    49564956        *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
    4957 
    49584957
    49594958        if(Jdet<0){
    … …  
    51495148}
    51505149/*}}}*/
    5151 /*FUNCTION Tria::GetParameterValue {{{1*/
     5150/*FUNCTION Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3) {{{1*/
    51525151void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){
    51535152       
    … …  
    51675166        /*Assign output pointers:*/
    51685167        *pp=p;
     5168}
     5169/*}}}*/
     5170/*FUNCTION Tria::GetParameterValue( {{{1*/
     5171void 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;
    51695227}
    51705228/*}}}*/

  • issm/trunk/src/c/objects/Elements/Tria.h

    r4402 r4422  
    150150                void      GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
    151151                void      GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
     152                void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
    152153                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    153154                void      GetSolutionFromInputsAdjoint(Vec solution);

  • issm/trunk/src/c/objects/FemModel.cpp

    r4402 r4422  
    143143
    144144/*Numerics: */
    145 /*FUNCTION FemModel::SetCurrentConfiguration{{{1*/
     145/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){{{1*/
    146146void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){
    147147
    … …  
    174174}
    175175/*}}}1*/
    176 /*FUNCTION FemModel::SetCurrentConfiguration{{{1*/
     176/*FUNCTION FemModel::SetCurrentConfiguration(int configuration_type){{{1*/
    177177void FemModel::SetCurrentConfiguration(int configuration_type){
    178178

  • issm/trunk/src/c/objects/Loads/Numericalflux.cpp

    r4405 r4422  
    432432
    433433                //Get vx, vy and v.n
    434                 trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum);
    435                 trias[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);
    436436               
    437437                UdotN=vx*normal[0]+vy*normal[1];
    … …  
    536536        GetNormal(&normal[0],xyz_list);
    537537
    538         /*Check wether it is an inflow or outflow BC*/
    539         trias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum);
    540         trias[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);
    541541        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    542542        if (UdotN<=0){
    … …  
    560560
    561561                //Get vx, vy and v.n
    562                 trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum);
    563                 trias[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);
    564564
    565565                UdotN=vx*normal[0]+vy*normal[1];
    … …  
    665665        GetNormal(&normal[0],xyz_list);
    666666
    667         /*Check wether it is an inflow or outflow BC*/
    668         trias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum);
    669         trias[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);
    670670        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    671671        if (UdotN>0){
    … …  
    689689
    690690                //Get vx, vy and v.n
    691                 trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum);
    692                 trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum);
    693                 trias[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);
    694694
    695695                UdotN=vx*normal[0]+vy*normal[1];
    … …  
    836836/*}}}*/
    837837/*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];
     838void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype){
    845839
    846840        /*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        }
    852866
    853867        /*Assign output pointers:*/
    854         *pp=p;
    855 }
    856 /*}}}*/
     868        *pvalue=value;
     869}
     870/*}}}*/

  • issm/trunk/src/c/objects/Loads/Numericalflux.h

    r4396 r4422  
    6868                void  GetDofList(int* doflist,int* pnumberofdofs);
    6969                void  GetNormal(double* normal,double xyz_list[4][3]);
    70                 void  GetParameterValue(double* pp, double* plist, double gauss_coord);
     70                void  GetParameterValue(double* pvalue, double gauss_coord,int enumtype);
    7171                void  CreateKMatrixInternal(Mat Kgg);
    7272                void  CreateKMatrixBoundary(Mat Kgg);
Note: See TracChangeset for help on using the changeset viewer.