Changeset 5808


Ignore:
Timestamp:
09/14/10 16:08:56 (15 years ago)
Author:
Mathieu Morlighem
Message:

New Icefront MUCH cleaner. Enjoy

Location:
issm/trunk/src/c
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r5806 r5808  
    43004300}
    43014301/*}}}*/
     4302/*FUNCTION Penta::GetZcoord {{{1*/
     4303double Penta::GetZcoord(GaussPenta* gauss){
     4304
     4305        int    i;
     4306        double xyz_list[NUMVERTICES][3];
     4307        double z_list[NUMVERTICES];
     4308        double z;
     4309
     4310        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4311        for(i=0;i<NUMVERTICES;i++){
     4312                z_list[i]=xyz_list[i][2];
     4313        }
     4314        PentaRef::GetParameterValue(&z,z_list,gauss);
     4315
     4316        return z;
     4317}
     4318/*}}}*/
    43024319/*FUNCTION Penta::GradjB {{{1*/
    43034320void  Penta::GradjB(Vec gradient){

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

    r5772 r5808  
    8080                int    GetNodeIndex(Node* node);
    8181                void   GetSolutionFromInputs(Vec solution);
     82                double GetZcoord(GaussPenta* gauss);
    8283                void   GetVectorFromInputs(Vec vector,int NameEnum);
    8384                void   Gradj(Vec gradient,int control_type);
    … …  
    146147                void      GetDofList(int** pdoflist,int approximation_enum,int setenum);
    147148                void      GetDofList1(int* doflist);
    148                 int       GetElementType(void);
    149                 void      GetParameterListOnVertices(double* pvalue,int enumtype);
    150                 void      GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
    151                 void      GetParameterValue(double* pvalue,Node* node,int enumtype);
     149                int     GetElementType(void);
     150                void    GetParameterListOnVertices(double* pvalue,int enumtype);
     151                void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
     152                void    GetParameterValue(double* pvalue,Node* node,int enumtype);
    152153                void      GetPhi(double* phi, double*  epsilon, double viscosity);
    153154                void      GetSolutionFromInputsDiagnosticHoriz(Vec solutiong);

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

    r5745 r5808  
    723723        z3=*(xyz_list+3*2+2);
    724724
     725        /*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
    725726        *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
    726727        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
    … …  
    935936        *(dl1dl6+NUMNODESP1*1+5)=1.0/SQRT3*(1.0+z)/2.0;
    936937        *(dl1dl6+NUMNODESP1*2+5)=0.5*A3;
     938}
     939/*}}}*/
     940/*FUNCTION PentaRef::GetQuadNodalFunctions {{{1*/
     941void PentaRef::GetQuadNodalFunctions(double* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4){
     942        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     943
     944        double BasisFunctions[6];
     945
     946        GetNodalFunctionsP1(&BasisFunctions[0],gauss);
     947
     948        ISSMASSERT(index1>=0 && index1<6);
     949        ISSMASSERT(index2>=0 && index2<6);
     950        ISSMASSERT(index3>=0 && index3<6);
     951        ISSMASSERT(index4>=0 && index4<6);
     952
     953        l1l4[0]=BasisFunctions[index1];
     954        l1l4[1]=BasisFunctions[index2];
     955        l1l4[2]=BasisFunctions[index3];
     956        l1l4[3]=BasisFunctions[index4];
     957
     958}
     959/*}}}*/
     960/*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{1*/
     961void PentaRef::GetQuadJacobianDeterminant(double* Jdet,double xyz_list[4][3],GaussPenta* gauss){
     962        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     963
     964        double x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
     965
     966        x1=xyz_list[0][0];
     967        y1=xyz_list[0][1];
     968        z1=xyz_list[0][2];
     969        x2=xyz_list[1][0];
     970        y2=xyz_list[1][1];
     971        z2=xyz_list[1][2];
     972        x3=xyz_list[2][0];
     973        y3=xyz_list[2][1];
     974        z3=xyz_list[2][2];
     975        x4=xyz_list[3][0];
     976        y4=xyz_list[3][1];
     977        z4=xyz_list[3][2];
     978
     979        /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/
     980        /*Area of a trabezoid = altitude * (base1 + base2)/2 */
     981        *Jdet= pow(pow(x2-x1,2.) + pow(y2-y1,2.),0.5) * (z4-z1 + z3-z2)/8;
     982        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     983
    937984}
    938985/*}}}*/

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

    r5743 r5808  
    2929                void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,GaussPenta* gauss);
    3030                void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,GaussPenta* gauss);
     31                void GetQuadNodalFunctions(double* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4);
     32                void GetQuadJacobianDeterminant(double*  Jdet, double xyz_list[4][3],GaussPenta* gauss);
    3133                void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss);
    3234                void GetJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);

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

    r5745 r5808  
    265265}
    266266/*}}}*/
    267 /*FUNCTION TriaRef::GetSegmentJacobian{{{1*/
    268 void TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss){
    269         /*The Jacobian is constant over the element, discard the gaussian points.
    270          * J is assumed to have been allocated.*/
    271 
     267/*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{1*/
     268void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){
     269        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     270         * J is assumed to have been allocated*/
    272271        double x1,y1,x2,y2;
    273272
    … …  
    277276        y2=*(xyz_list+3*1+1);
    278277
    279         *J=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.));
    280 }
    281 /*}}}*/
    282 /*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{1*/
    283 void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){
    284         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    285          * J is assumed to have been allocated*/
    286 
    287         GetSegmentJacobian(Jdet,xyz_list,gauss);
     278        *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.));
    288279        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
    289280

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

    r5743 r5808  
    3131                void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof);
    3232                void GetJacobian(double* J, double* xyz_list,GaussTria* gauss);
    33                 void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss);
    3433                void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss);
    3534                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss);

  • issm/trunk/src/c/objects/Gauss/GaussPenta.cpp

    r5797 r5808  
    179179
    180180        /*Quads: get the gauss points using the product of two line rules  */
    181         if(index1==0 && index2==1 && index3==3 && index4==4){
     181        if(index1==0 && index2==1 && index3==4 && index4==3){
    182182                for(i=0;i<order_horiz;i++){
    183183                        for(j=0;j<order_vert;j++){
    184                                 coords1[i*order_horiz+j]=  0.5*(1-seg_horiz_coords[i]);
    185                                 coords2[i*order_horiz+j]=1-0.5*(1-seg_horiz_coords[i]);
    186                                 coords3[i*order_horiz+j]=0.0;
    187                                 coords4[i*order_horiz+j]=seg_vert_coords[j];
    188                                 weights[i*order_horiz+j]=seg_horiz_weights[i]*seg_vert_weights[j];
     184                                coords1[i*order_vert+j]=  0.5*(1-seg_horiz_coords[i]);
     185                                coords2[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
     186                                coords3[i*order_vert+j]=0.0;
     187                                coords4[i*order_vert+j]=seg_vert_coords[j];
     188                                weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
    189189                        }
    190190                }
    191191        }
    192         else if(index1==1 && index2==2 && index3==4 && index4==5){
     192        else if(index1==1 && index2==2 && index3==5 && index4==4){
    193193                for(i=0;i<order_horiz;i++){
    194194                        for(j=0;j<order_vert;j++){
    195                                 coords1[i*order_horiz+j]=0.0;
    196                                 coords2[i*order_horiz+j]=  0.5*(1-seg_horiz_coords[i]);
    197                                 coords3[i*order_horiz+j]=1-0.5*(1-seg_horiz_coords[i]);
    198                                 coords4[i*order_horiz+j]=seg_vert_coords[j];
    199                                 weights[i*order_horiz+j]=seg_horiz_weights[i]*seg_vert_weights[j];
     195                                coords1[i*order_vert+j]=0.0;
     196                                coords2[i*order_vert+j]=  0.5*(1-seg_horiz_coords[i]);
     197                                coords3[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
     198                                coords4[i*order_vert+j]=seg_vert_coords[j];
     199                                weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
    200200                        }
    201201                }
    202202        }
    203         else if(index1==2 && index2==0 && index3==5 && index4==3){
     203        else if(index1==2 && index2==0 && index3==3 && index4==5){
    204204                for(i=0;i<order_horiz;i++){
    205205                        for(j=0;j<order_vert;j++){
    206                                 coords1[i*order_horiz+j]=1-0.5*(1-seg_horiz_coords[i]);
    207                                 coords2[i*order_horiz+j]=0.0;
    208                                 coords3[i*order_horiz+j]=  0.5*(1-seg_horiz_coords[i]);
    209                                 coords4[i*order_horiz+j]=seg_vert_coords[j];
    210                                 weights[i*order_horiz+j]=seg_horiz_weights[i]*seg_vert_weights[j];
     206                                coords1[i*order_vert+j]=1-0.5*(1-seg_horiz_coords[i]);
     207                                coords2[i*order_vert+j]=0.0;
     208                                coords3[i*order_vert+j]=  0.5*(1-seg_horiz_coords[i]);
     209                                coords4[i*order_vert+j]=seg_vert_coords[j];
     210                                weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
    211211                        }
    212212                }
    … …  
    216216        }
    217217
     218        /*clean-up*/
     219        xfree((void**)&seg_horiz_coords);
     220        xfree((void**)&seg_horiz_weights);
     221        xfree((void**)&seg_vert_coords);
     222        xfree((void**)&seg_vert_weights);
    218223}
    219224/*}}}*/

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

    r5772 r5808  
    457457        gravity  =matpar->GetG();
    458458
    459         GetNormal(&normal[0],xyz_list);
     459        GetSegmentNormal(&normal[0],xyz_list);
    460460
    461461        index1=tria->GetNodeIndex(nodes[0]);
    … …  
    542542void Icefront::CreatePVectorDiagnosticPattyn( Vec pg){
    543543
    544         int i,j;
    545 
    546         const int numnodes=NUMVERTICESQUA;
    547         const int numdofs=numnodes*NDOF2;
    548         int*      doflist=NULL;
    549 
    550         double xyz_list[numnodes][3];
    551         double xyz_list_quad[numnodes+1][3]; //center node
    552 
    553         double thickness_list[6]; //from penta element
    554         double thickness_list_quad[6]; //selection of 4 thickness from penta elements
    555 
    556         double bed_list[6]; //from penta element
    557         double bed_list_quad[6]; //selection of 4 beds from penta elements
    558 
    559         /*Objects: */
    560         double  pe_g[numdofs]={0.0};
    561         int    element_type;
    562                
    563         /*quad grids: */
    564         int grid1,grid2,grid3,grid4;
    565 
    566         /*normals: */
    567         double normal1[3];
    568         double normal2[3];
    569         double normal3[3];
    570         double normal4[3];
    571         double normal_norm;
    572         double v15[3];
    573         double v25[3];
    574         double v35[3];
    575         double v45[3];
    576 
    577 
    578         element_type=element->Enum();
    579 
    580         /*check icefront is associated to a pentaelem: */
    581         if(element_type!=PentaEnum){
    582                 ISSMERROR(" quad icefront is associated to a TriaElem!");
    583         }
     544        /*Constants*/
     545        const int numdofs = NUMVERTICESQUA *NDOF2;
     546
     547        /*Intermediaries*/
     548        int         i,j,ig,index1,index2,index3,index4;
     549        int         fill;
     550        double      surface,pressure,ice_pressure,rho_water,rho_ice,gravity;
     551        double      water_pressure,air_pressure;
     552        double      Jdet,z_g;
     553        double      xyz_list[NUMVERTICESQUA][3];
     554        double      normal[3];
     555        double      pe_g[numdofs] = {0.0};
     556        double      l1l4[4];
     557        int        *doflist = NULL;
     558        GaussPenta *gauss = NULL;
     559        Penta      *penta = NULL;
    584560
    585561        /* Get dof list and node coordinates: */
     562        penta=(Penta*)element;
    586563        GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    587         GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    588        
    589         //Identify which grids are comprised in the quad:
    590         grid1=element->GetNodeIndex(nodes[0]);
    591         grid2=element->GetNodeIndex(nodes[1]);
    592         grid3=element->GetNodeIndex(nodes[2]);
    593         grid4=element->GetNodeIndex(nodes[3]);
    594 
    595         /*Build new xyz, bed and thickness lists for grids 1 to 4: */
    596         for(i=0;i<4;i++){
    597                 for(j=0;j<3;j++){
    598                         xyz_list_quad[i][j]=xyz_list[i][j];
    599                 }
    600         }
    601 
    602         element->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum);
    603         element->GetParameterListOnVertices(&bed_list[0],BedEnum);
    604 
    605         thickness_list_quad[0]=thickness_list[grid1];
    606         thickness_list_quad[1]=thickness_list[grid2];
    607         thickness_list_quad[2]=thickness_list[grid3];
    608         thickness_list_quad[3]=thickness_list[grid4];
    609 
    610         bed_list_quad[0]=bed_list[grid1];
    611         bed_list_quad[1]=bed_list[grid2];
    612         bed_list_quad[2]=bed_list[grid3];
    613         bed_list_quad[3]=bed_list[grid4];
    614 
    615         //Create a new grid in the midle of the quad and add it at the end of the list
    616         xyz_list_quad[4][0] = (xyz_list_quad[0][0]+xyz_list_quad[1][0]+xyz_list_quad[2][0]+xyz_list_quad[3][0])/4.0;
    617         xyz_list_quad[4][1] = (xyz_list_quad[0][1]+xyz_list_quad[1][1]+xyz_list_quad[2][1]+xyz_list_quad[3][1])/4.0;
    618         xyz_list_quad[4][2] = (xyz_list_quad[0][2]+xyz_list_quad[1][2]+xyz_list_quad[2][2]+xyz_list_quad[3][2])/4.0;
    619 
    620         //Compute four normals (the quad is divided into four triangles)
    621         v15[0]=xyz_list_quad[0][0]-xyz_list_quad[4][0];
    622         v15[1]=xyz_list_quad[0][1]-xyz_list_quad[4][1];
    623         v15[2]=xyz_list_quad[0][2]-xyz_list_quad[4][2];
    624        
    625         v25[0]=xyz_list_quad[1][0]-xyz_list_quad[4][0];
    626         v25[1]=xyz_list_quad[1][1]-xyz_list_quad[4][1];
    627         v25[2]=xyz_list_quad[1][2]-xyz_list_quad[4][2];
    628        
    629         v35[0]=xyz_list_quad[2][0]-xyz_list_quad[4][0];
    630         v35[1]=xyz_list_quad[2][1]-xyz_list_quad[4][1];
    631         v35[2]=xyz_list_quad[2][2]-xyz_list_quad[4][2];
    632        
    633         v45[0]=xyz_list_quad[3][0]-xyz_list_quad[4][0];
    634         v45[1]=xyz_list_quad[3][1]-xyz_list_quad[4][1];
    635         v45[2]=xyz_list_quad[3][2]-xyz_list_quad[4][2];
    636 
    637         cross(normal1,v15,v25);
    638         cross(normal2,v25,v35);
    639         cross(normal3,v35,v45);
    640         cross(normal4,v45,v15);
    641 
    642         normal_norm=norm(normal1);normal1[0]=normal1[0]/normal_norm;normal1[1]=normal1[1]/normal_norm;normal1[2]=normal1[2]/normal_norm;
    643         normal_norm=norm(normal2);normal2[0]=normal2[0]/normal_norm;normal2[1]=normal2[1]/normal_norm;normal2[2]=normal2[2]/normal_norm;
    644         normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm;
    645         normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm;
    646 
    647         //Compute load contribution for this quad:
    648         QuadPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0]);
    649 
    650         /*Plug pe_g into vector: */
     564        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA);
     565
     566        /*Recover inputs and parameters */
     567        Input* surface_input  =penta->inputs->GetInput(SurfaceEnum);   ISSMASSERT(surface_input);
     568        inputs->GetParameterValue(&fill,FillEnum);
     569        rho_water=matpar->GetRhoWater();
     570        rho_ice  =matpar->GetRhoIce();
     571        gravity  =matpar->GetG();
     572        GetQuadNormal(&normal[0],xyz_list);
     573
     574        /*Identify which nodes are in the quad: */
     575        index1=element->GetNodeIndex(nodes[0]);
     576        index2=element->GetNodeIndex(nodes[1]);
     577        index3=element->GetNodeIndex(nodes[2]);
     578        index4=element->GetNodeIndex(nodes[3]);
     579
     580        /* Start  looping on the number of gaussian points: */
     581        double zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
     582        double zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
     583        if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity
     584        else                 gauss=new GaussPenta(index1,index2,index3,index4,3,3);
     585        for(ig=gauss->begin();ig<gauss->end();ig++){
     586
     587                gauss->GaussPoint(ig);
     588
     589                penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
     590                penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
     591                z_g=penta->GetZcoord(gauss);
     592                surface_input->GetParameterValue(&surface,gauss);
     593
     594                switch(fill){
     595                        case WaterEnum:
     596                                water_pressure=rho_water*gravity*min(0,z_g);//0 if the gaussian point is above water level
     597                                break;
     598                        case AirEnum:
     599                                water_pressure=0;
     600                                break;
     601                        default:
     602                                ISSMERROR("fill type %s not supported yet",EnumToString(fill));
     603                }
     604                ice_pressure=rho_ice*gravity*(surface-z_g);
     605                air_pressure=0;
     606                pressure = ice_pressure + water_pressure + air_pressure;
     607
     608                for(i=0;i<NUMVERTICESQUA;i++){
     609                        for(j=0;j<NDOF2;j++){
     610                                pe_g[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
     611                        }
     612                }
     613        }
     614
    651615        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    652616
    653617        /*Free ressources:*/
     618        delete gauss;
    654619        xfree((void**)&doflist);
    655 
    656 }
    657 /*}}}*/
    658 /*FUNCTION Icefront::CreatePVectorDiagnosticStokes {{{1*/
     620}
     621/*}}}*/
     622/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{1*/
    659623void Icefront::CreatePVectorDiagnosticStokes( Vec pg){
    660624
    661         int i,j;
    662 
    663         const int numnodes=4;
    664         const int numdofs=numnodes*NDOF4;
    665         int*      doflist=NULL;
    666 
    667         double xyz_list[numnodes][3];
    668         double xyz_list_quad[numnodes+1][3]; //center node
    669 
    670         double thickness_list[6]; //from penta element
    671         double thickness_list_quad[6]; //selection of 4 thickness from penta elements
    672 
    673         double bed_list[6]; //from penta element
    674         double bed_list_quad[6]; //selection of 4 beds from penta elements
    675 
    676         /*Objects: */
    677         double  pe_g[numdofs]={0.0};
    678         Penta*   penta=NULL;
    679 
    680         /*quad grids: */
    681         int grid1,grid2,grid3,grid4;
    682 
    683         /*normals: */
    684         double normal1[3];
    685         double normal2[3];
    686         double normal3[3];
    687         double normal4[3];
    688         double normal_norm;
    689         double v15[3];
    690         double v25[3];
    691         double v35[3];
    692         double v45[3];
    693         int approximation;
    694 
    695         /*check icefront is associated to a pentaelem: */
    696         if(element->Enum()!=PentaEnum) ISSMERROR("Only Penta supported yet");
     625        /*Constants*/
     626        const int numdofs = NUMVERTICESQUA *NDOF4;
     627
     628        /*Intermediaries*/
     629        int         i,j,ig,index1,index2,index3,index4;
     630        int         fill;
     631        double      pressure,rho_water,gravity;
     632        double      water_pressure,air_pressure;
     633        double      Jdet,z_g;
     634        double      xyz_list[NUMVERTICESQUA][3];
     635        double      normal[3];
     636        double      pe_g[numdofs] = {0.0};
     637        double      l1l4[4];
     638        int        *doflist = NULL;
     639        GaussPenta *gauss = NULL;
     640        Penta      *penta = NULL;
     641
     642        /* Get dof list and node coordinates: */
    697643        penta=(Penta*)element;
    698 
    699         /*If not Stokes, return*/
    700         penta->inputs->GetParameterValue(&approximation,ApproximationEnum);
    701         if(approximation!=StokesApproximationEnum) return;
    702 
    703         /* Get dof list and node coordinates: */
    704644        GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    705         GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
    706        
    707         //Identify which grids are comprised in the quad:
    708         grid1=element->GetNodeIndex(nodes[0]);
    709         grid2=element->GetNodeIndex(nodes[1]);
    710         grid3=element->GetNodeIndex(nodes[2]);
    711         grid4=element->GetNodeIndex(nodes[3]);
    712 
    713         /*Build new xyz, bed and thickness lists for grids 1 to 4: */
    714         for(i=0;i<4;i++){
    715                 for(j=0;j<3;j++){
    716                         xyz_list_quad[i][j]=xyz_list[i][j];
    717                 }
    718         }
    719 
    720         element->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum);
    721         element->GetParameterListOnVertices(&bed_list[0],BedEnum);
    722 
    723         thickness_list_quad[0]=thickness_list[grid1];
    724         thickness_list_quad[1]=thickness_list[grid2];
    725         thickness_list_quad[2]=thickness_list[grid3];
    726         thickness_list_quad[3]=thickness_list[grid4];
    727 
    728         bed_list_quad[0]=bed_list[grid1];
    729         bed_list_quad[1]=bed_list[grid2];
    730         bed_list_quad[2]=bed_list[grid3];
    731         bed_list_quad[3]=bed_list[grid4];
    732 
    733         //Create a new grid in the midle of the quad and add it at the end of the list
    734         xyz_list_quad[4][0] = (xyz_list_quad[0][0]+xyz_list_quad[1][0]+xyz_list_quad[2][0]+xyz_list_quad[3][0])/4.0;
    735         xyz_list_quad[4][1] = (xyz_list_quad[0][1]+xyz_list_quad[1][1]+xyz_list_quad[2][1]+xyz_list_quad[3][1])/4.0;
    736         xyz_list_quad[4][2] = (xyz_list_quad[0][2]+xyz_list_quad[1][2]+xyz_list_quad[2][2]+xyz_list_quad[3][2])/4.0;
    737 
    738         //Compute four normals (the quad is divided into four triangles)
    739         v15[0]=xyz_list_quad[0][0]-xyz_list_quad[4][0];
    740         v15[1]=xyz_list_quad[0][1]-xyz_list_quad[4][1];
    741         v15[2]=xyz_list_quad[0][2]-xyz_list_quad[4][2];
    742        
    743         v25[0]=xyz_list_quad[1][0]-xyz_list_quad[4][0];
    744         v25[1]=xyz_list_quad[1][1]-xyz_list_quad[4][1];
    745         v25[2]=xyz_list_quad[1][2]-xyz_list_quad[4][2];
    746        
    747         v35[0]=xyz_list_quad[2][0]-xyz_list_quad[4][0];
    748         v35[1]=xyz_list_quad[2][1]-xyz_list_quad[4][1];
    749         v35[2]=xyz_list_quad[2][2]-xyz_list_quad[4][2];
    750        
    751         v45[0]=xyz_list_quad[3][0]-xyz_list_quad[4][0];
    752         v45[1]=xyz_list_quad[3][1]-xyz_list_quad[4][1];
    753         v45[2]=xyz_list_quad[3][2]-xyz_list_quad[4][2];
    754 
    755         cross(normal1,v15,v25);
    756         cross(normal2,v25,v35);
    757         cross(normal3,v35,v45);
    758         cross(normal4,v45,v15);
    759 
    760         normal_norm=norm(normal1);normal1[0]=normal1[0]/normal_norm;normal1[1]=normal1[1]/normal_norm;normal1[2]=normal1[2]/normal_norm;
    761         normal_norm=norm(normal2);normal2[0]=normal2[0]/normal_norm;normal2[1]=normal2[1]/normal_norm;normal2[2]=normal2[2]/normal_norm;
    762         normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm;
    763         normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm;
    764 
    765 
    766         //Compute load contribution for this quad:
    767         QuadPressureLoadStokes(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0]);
    768 
    769         /*Plug pe_g into vector: */
     645        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA);
     646       
     647        /*Recover inputs and parameters */
     648        inputs->GetParameterValue(&fill,FillEnum);
     649        rho_water=matpar->GetRhoWater();
     650        gravity  =matpar->GetG();
     651        GetQuadNormal(&normal[0],xyz_list);
     652
     653        /*Identify which nodes are in the quad: */
     654        index1=element->GetNodeIndex(nodes[0]);
     655        index2=element->GetNodeIndex(nodes[1]);
     656        index3=element->GetNodeIndex(nodes[2]);
     657        index4=element->GetNodeIndex(nodes[3]);
     658
     659        /* Start  looping on the number of gaussian points: */
     660        double zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
     661        double zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
     662        if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity
     663        else                 gauss=new GaussPenta(index1,index2,index3,index4,3,3);
     664        for(ig=gauss->begin();ig<gauss->end();ig++){
     665
     666                gauss->GaussPoint(ig);
     667
     668                penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
     669                penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
     670                z_g=penta->GetZcoord(gauss);
     671
     672                switch(fill){
     673                        case WaterEnum:
     674                                water_pressure=rho_water*gravity*min(0,z_g);//0 if the gaussian point is above water level
     675                                break;
     676                        case AirEnum:
     677                                water_pressure=0;
     678                                break;
     679                        default:
     680                                ISSMERROR("fill type %s not supported yet",EnumToString(fill));
     681                }
     682                air_pressure=0;
     683                pressure = water_pressure + air_pressure; //no ice pressure fore Stokes
     684
     685                for(i=0;i<NUMVERTICESQUA;i++){
     686                        for(j=0;j<NDOF4;j++){
     687                                if(j<3)  pe_g[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
     688                                else     pe_g[i*NDOF4+j]+=0; //pressure term
     689                        }
     690                }
     691        }
     692
    770693        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    771694
    772695        /*Free ressources:*/
     696        delete gauss;
    773697        xfree((void**)&doflist);
    774 
    775698}
    776699/*}}}*/
    … …  
    949872}
    950873/*}}}*/
    951 /*FUNCTION Icefront::QuadPressureLoad {{{1*/
    952 void Icefront::QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    953                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list){
    954        
    955        
    956         //The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
    957         //
    958         //   grid 4 +-----------------+ grid 3
    959         //          |\2            1 /|
    960         //          |1\    tria3    /2|
    961         //          |  \           /  |
    962         //          |   \         /   |
    963         //          |    \       /    |
    964         //          |     \     /     |
    965         //          |      \ 3 /      |
    966         //          |tria4  \ / 3     |
    967         //          |      3 \grid5   |
    968         //          |       / \       |
    969         //          |      / 3 \ tria2|
    970         //          |     /     \     |
    971         //          |    /       \    |
    972         //          |   /         \   |
    973         //          |  /   tria1   \  |
    974         //          |2/1           2\1|
    975         //   grid1  +-----------------+ grid 2
    976         //
    977         //
    978 
    979         int      i,j;
    980 
    981         double nx[4];
    982         double ny[4];
    983 
    984         /* gaussian points: */
    985         int     ig;
    986         GaussTria *gauss=NULL;
    987 
    988         double  surface_list[5];
    989         double  x[5];
    990         double  y[5];
    991         double  z[5];
    992         double l12,l14,l15,l23,l25,l34,l35,l45;
    993         double cos_theta_triangle1,cos_theta_triangle2,cos_theta_triangle3,cos_theta_triangle4;
    994 
    995         double xyz_tria[12][2];
    996         double l1l2l3_tria[3];
    997         double r_tria,s_tria;
    998         double r_quad[4];
    999         double s_quad[4];
    1000         double l1l4_tria[4][4];
    1001 
    1002         double J[4];
    1003         double z_g[4];
    1004         double ice_pressure_tria[4];
    1005         double air_pressure_tria;
    1006         double water_level_above_g_tria;
    1007         double water_pressure_tria;
    1008         double pressure_tria[4];
    1009         int    fill;
    1010 
    1011 
    1012         /*To use tria functionalities: */
    1013         Tria* tria=NULL;
    1014        
    1015         /*Recover inputs: */
    1016         inputs->GetParameterValue(&fill,FillEnum);
    1017 
    1018         /*Initialize fake tria: */
    1019         tria=new Tria();
    1020 
    1021         //Build the four normal vectors
    1022         nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0];
    1023         ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
    1024        
    1025         //Recover the surface of the four nodes
    1026         for(i=0;i<4;i++){
    1027                 surface_list[i]=thickness_list[i]+bed_list[i];
    1028         }
    1029         //Add surface sor the fifth point (average of the surfaces)
    1030         surface_list[4]=(surface_list[0]+surface_list[1]+surface_list[2]+surface_list[3])/4.0;
    1031 
    1032         //Recover grid coordinates
    1033         for(i=0;i<5;i++){
    1034                 x[i]=*(xyz_list+3*i+0);
    1035                 y[i]=*(xyz_list+3*i+1);
    1036                 z[i]=*(xyz_list+3*i+2);
    1037         }
    1038        
    1039         //Build triangles in a 2D plan before using reference elements
    1040         l12=pow((x[1]-x[0]),2)+pow((y[1]-y[0]),2)+pow((z[1]-z[0]),2);
    1041         l14=pow((x[3]-x[0]),2)+pow((y[3]-y[0]),2)+pow((z[3]-z[0]),2);
    1042         l15=pow((x[4]-x[0]),2)+pow((y[4]-y[0]),2)+pow((z[4]-z[0]),2);
    1043         l23=pow((x[2]-x[1]),2)+pow((y[2]-y[1]),2)+pow((z[2]-z[1]),2);
    1044         l25=pow((x[4]-x[1]),2)+pow((y[4]-y[1]),2)+pow((z[4]-z[1]),2);
    1045         l34=pow((x[3]-x[2]),2)+pow((y[3]-y[2]),2)+pow((z[3]-z[2]),2);
    1046         l35=pow((x[4]-x[2]),2)+pow((y[4]-y[2]),2)+pow((z[4]-z[2]),2);
    1047         l45=pow((x[4]-x[3]),2)+pow((y[4]-y[3]),2)+pow((z[4]-z[3]),2);
    1048 
    1049         cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
    1050         cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
    1051         cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
    1052         cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
    1053 
    1054         //First triangle : nodes 1, 2 and 5
    1055         xyz_tria[0][0]=0;
    1056         xyz_tria[0][1]=0;
    1057         xyz_tria[1][0]=sqrt(l12);
    1058         xyz_tria[1][1]=0;
    1059         xyz_tria[2][0]=cos_theta_triangle1*sqrt(l15);
    1060         xyz_tria[2][1]=sqrt(1-pow(cos_theta_triangle1,2))*sqrt(l15);
    1061 
    1062         //Second triangle : nodes 2, 3 and 5
    1063         xyz_tria[3][0]=0;
    1064         xyz_tria[3][1]=0;
    1065         xyz_tria[4][0]=sqrt(l23);
    1066         xyz_tria[4][1]=0;
    1067         xyz_tria[5][0]=cos_theta_triangle2*sqrt(l25);
    1068         xyz_tria[5][1]=sqrt(1-pow(cos_theta_triangle2,2))*sqrt(l25);
    1069        
    1070         //Third triangle : nodes 3, 4 and 5
    1071         xyz_tria[6][0]=0;
    1072         xyz_tria[6][1]=0;
    1073         xyz_tria[7][0]=sqrt(l34);
    1074         xyz_tria[7][1]=0;
    1075         xyz_tria[8][0]=cos_theta_triangle3*sqrt(l35);
    1076         xyz_tria[8][1]=sqrt(1-pow(cos_theta_triangle3,2))*sqrt(l35);
    1077 
    1078         //Fourth triangle : nodes 4, 1 and 5
    1079         xyz_tria[9][0]=0;
    1080         xyz_tria[9][1]=0;
    1081         xyz_tria[10][0]=sqrt(l14);
    1082         xyz_tria[10][1]=0;
    1083         xyz_tria[11][0]=cos_theta_triangle4*sqrt(l45);
    1084         xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
    1085 
    1086         /* Start  looping on the number of gaussian points: */
    1087         gauss=new GaussTria(2);
    1088         for(ig=gauss->begin();ig<gauss->end();ig++){
    1089 
    1090                 gauss->GaussPoint(ig);
    1091 
    1092                 tria->GetNodalFunctions(l1l2l3_tria,gauss);
    1093 
    1094                 //Get the coordinates of gauss point for each triangle in penta/quad
    1095                 r_tria=gauss->coord2-gauss->coord1;
    1096                 s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3);
    1097 
    1098                 //Coordinates of gauss points in the reference triangle
    1099                 r_quad[0]=r_tria;
    1100                 s_quad[0]=1/sqrt((double)3)*s_tria-2/3;
    1101                 r_quad[1]=-1/sqrt((double)3)*s_tria+2/3;
    1102                 s_quad[1]=r_tria;
    1103                 r_quad[2]=-r_tria;
    1104                 s_quad[2]=-1/sqrt((double)3)*s_tria+2/3;
    1105                 r_quad[3]=1/sqrt((double)3)*s_tria-2/3;
    1106                 s_quad[3]=-r_tria;
    1107 
    1108                 //Get the nodal function of the quad for the gauss points of each triangle
    1109                 for(i=0;i<4;i++){
    1110                         l1l4_tria[i][0]=1.0/4.0*(
    1111                                         (r_quad[i]-1.0)*(s_quad[i]-1.0)
    1112                                         );
    1113                        
    1114                         l1l4_tria[i][1]=1.0/4.0*(
    1115                                          -(r_quad[i]+1.0)*(s_quad[i]-1.0)
    1116                                         );
    1117 
    1118                         l1l4_tria[i][2]=1.0/4.0*(
    1119                                         (r_quad[i]+1.0)*(s_quad[i]+1.0)
    1120                                         );
    1121                        
    1122                         l1l4_tria[i][3]=1.0/4.0*(
    1123                                          -(r_quad[i]-1.0)*(s_quad[i]+1.0)
    1124                                         );
    1125 
    1126                 } //for(i=0;i<4;i++)
    1127                
    1128                
    1129                 //Compute jacobian of each triangle
    1130                 for (i=0;i<4;i++){
    1131                         double complete_list[3][3]; //a third coordinate is needed for the jacobian determinant calculation, here it is zero
    1132                         for(j=0;j<3;j++){
    1133                                 complete_list[j][0]=xyz_tria[3*i+j][0];
    1134                                 complete_list[j][1]=xyz_tria[3*i+j][1];
    1135                                 complete_list[j][2]=0;
    1136                         }
    1137                         tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss);
    1138                 }
    1139 
    1140                 //Calculation of the z coordinate for the gaussian point ig for each triangle
    1141                 z_g[0]=z[0]*l1l2l3_tria[0]+z[1]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1142                 z_g[1]=z[1]*l1l2l3_tria[0]+z[2]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1143                 z_g[2]=z[2]*l1l2l3_tria[0]+z[3]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1144                 z_g[3]=z[3]*l1l2l3_tria[0]+z[0]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1145                
    1146                 //Loop on the triangles
    1147                 for(i=0;i<4;i++){
    1148 
    1149                         //Loop on the grids of the quad
    1150                         //Calculate the ice pressure
    1151                         for(j=0;j<4;j++){
    1152                                 ice_pressure_tria[j]=gravity*rho_ice*(surface_list[j]-z_g[i]);
    1153                         }
    1154                         air_pressure_tria=0;
    1155 
    1156                         //Now deal with water pressure:
    1157                         if(fill==WaterEnum){ //icefront ends in water
    1158                                 water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level
    1159                                 water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
    1160                         }
    1161                         else if(fill==AirEnum){
    1162                                 water_pressure_tria=0;
    1163                         }
    1164                         else{
    1165                                 ISSMERROR("QuadPressureLoad error message: unknow fill type for icefront boundary condition");
    1166                         }
    1167 
    1168                         //Add pressure from triangle i
    1169                         //Loop on the grids of the quad
    1170                         for(j=0;j<4;j++){
    1171                                 pressure_tria[j] = ice_pressure_tria[j] + water_pressure_tria + air_pressure_tria;
    1172                         }
    1173 
    1174 
    1175                         pe_g[0]+= J[i]*gauss->weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i];
    1176                         pe_g[1]+= J[i]*gauss->weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i];
    1177                         pe_g[2]+= J[i]*gauss->weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i];
    1178                         pe_g[3]+= J[i]*gauss->weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i];
    1179                         pe_g[4]+= J[i]*gauss->weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i];
    1180                         pe_g[5]+= J[i]*gauss->weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i];
    1181                         pe_g[6]+= J[i]*gauss->weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i];
    1182                         pe_g[7]+= J[i]*gauss->weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i];
    1183 
    1184                        
    1185 
    1186                 } //for(i=0;i<4;i++)
    1187         }
    1188        
    1189         delete tria;
    1190         delete gauss;
    1191 }
    1192 /*}}}*/
    1193 /*FUNCTION Icefront::QuadPressureLoadStokes {{{1*/
    1194 void Icefront::QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    1195                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list){
    1196        
    1197        
    1198         //The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
    1199         //
    1200         //   grid 4 +-----------------+ grid 3
    1201         //          |\2            1 /|
    1202         //          |1\    tria3    /2|
    1203         //          |  \           /  |
    1204         //          |   \         /   |
    1205         //          |    \       /    |
    1206         //          |     \     /     |
    1207         //          |      \ 3 /      |
    1208         //          |tria4  \ / 3     |
    1209         //          |      3 \grid5   |
    1210         //          |       / \       |
    1211         //          |      / 3 \ tria2|
    1212         //          |     /     \     |
    1213         //          |    /       \    |
    1214         //          |   /         \   |
    1215         //          |  /   tria1   \  |
    1216         //          |2/1           2\1|
    1217         //   grid1  +-----------------+ grid 2
    1218         //
    1219         //
    1220 
    1221         int      i,j;
    1222 
    1223         double nx[4];
    1224         double ny[4];
    1225         double nz[4];
    1226 
    1227         /* gaussian points: */
    1228         int     ig;
    1229         GaussTria* gauss=NULL;
    1230 
    1231         double  surface_list[5];
    1232         double  x[5];
    1233         double  y[5];
    1234         double  z[5];
    1235         double l12,l14,l15,l23,l25,l34,l35,l45;
    1236         double cos_theta_triangle1,cos_theta_triangle2,cos_theta_triangle3,cos_theta_triangle4;
    1237 
    1238         double xyz_tria[12][2];
    1239         double l1l2l3_tria[3];
    1240         double r_tria,s_tria;
    1241         double r_quad[4];
    1242         double s_quad[4];
    1243         double l1l4_tria[4][4];
    1244 
    1245         double J[4];
    1246         double z_g[4];
    1247         double air_pressure_tria;
    1248         double water_level_above_g_tria;
    1249         double water_pressure_tria;
    1250         double pressure_tria;
    1251         int    fill;
    1252 
    1253         /*To use tria functionalities: */
    1254         Tria* tria=NULL;
    1255 
    1256         /*Recover inputs: */
    1257         inputs->GetParameterValue(&fill,FillEnum);
    1258 
    1259         /*Initialize fake tria: */
    1260         tria=new Tria();
    1261 
    1262         //Build the four normal vectors
    1263         nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0];
    1264         ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1];
    1265         nz[0]=normal1[2]; nz[1]=normal2[2]; nz[2]=normal3[2]; nz[3]=normal4[2];
    1266        
    1267         //Recover the surface of the four nodes
    1268         for(i=0;i<4;i++){
    1269                 surface_list[i]=thickness_list[i]+bed_list[i];
    1270         }
    1271         //Add surface sor the fifth point (average of the surfaces)
    1272         surface_list[4]=(surface_list[0]+surface_list[1]+surface_list[2]+surface_list[3])/4.0;
    1273 
    1274         //Recover grid coordinates
    1275         for(i=0;i<5;i++){
    1276                 x[i]=*(xyz_list+3*i+0);
    1277                 y[i]=*(xyz_list+3*i+1);
    1278                 z[i]=*(xyz_list+3*i+2);
    1279         }
    1280        
    1281         //Build triangles in a 2D plan before using reference elements
    1282         l12=pow((x[1]-x[0]),2)+pow((y[1]-y[0]),2)+pow((z[1]-z[0]),2);
    1283         l14=pow((x[3]-x[0]),2)+pow((y[3]-y[0]),2)+pow((z[3]-z[0]),2);
    1284         l15=pow((x[4]-x[0]),2)+pow((y[4]-y[0]),2)+pow((z[4]-z[0]),2);
    1285         l23=pow((x[2]-x[1]),2)+pow((y[2]-y[1]),2)+pow((z[2]-z[1]),2);
    1286         l25=pow((x[4]-x[1]),2)+pow((y[4]-y[1]),2)+pow((z[4]-z[1]),2);
    1287         l34=pow((x[3]-x[2]),2)+pow((y[3]-y[2]),2)+pow((z[3]-z[2]),2);
    1288         l35=pow((x[4]-x[2]),2)+pow((y[4]-y[2]),2)+pow((z[4]-z[2]),2);
    1289         l45=pow((x[4]-x[3]),2)+pow((y[4]-y[3]),2)+pow((z[4]-z[3]),2);
    1290 
    1291         cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
    1292         cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
    1293         cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
    1294         cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));
    1295 
    1296         //First triangle : nodes 1, 2 and 5
    1297         xyz_tria[0][0]=0;
    1298         xyz_tria[0][1]=0;
    1299         xyz_tria[1][0]=sqrt(l12);
    1300         xyz_tria[1][1]=0;
    1301         xyz_tria[2][0]=cos_theta_triangle1*sqrt(l15);
    1302         xyz_tria[2][1]=sqrt(1-pow(cos_theta_triangle1,2))*sqrt(l15);
    1303 
    1304         //Second triangle : nodes 2, 3 and 5
    1305         xyz_tria[3][0]=0;
    1306         xyz_tria[3][1]=0;
    1307         xyz_tria[4][0]=sqrt(l23);
    1308         xyz_tria[4][1]=0;
    1309         xyz_tria[5][0]=cos_theta_triangle2*sqrt(l25);
    1310         xyz_tria[5][1]=sqrt(1-pow(cos_theta_triangle2,2))*sqrt(l25);
    1311        
    1312         //Third triangle : nodes 3, 4 and 5
    1313         xyz_tria[6][0]=0;
    1314         xyz_tria[6][1]=0;
    1315         xyz_tria[7][0]=sqrt(l34);
    1316         xyz_tria[7][1]=0;
    1317         xyz_tria[8][0]=cos_theta_triangle3*sqrt(l35);
    1318         xyz_tria[8][1]=sqrt(1-pow(cos_theta_triangle3,2))*sqrt(l35);
    1319 
    1320         //Fourth triangle : nodes 4, 1 and 5
    1321         xyz_tria[9][0]=0;
    1322         xyz_tria[9][1]=0;
    1323         xyz_tria[10][0]=sqrt(l14);
    1324         xyz_tria[10][1]=0;
    1325         xyz_tria[11][0]=cos_theta_triangle4*sqrt(l45);
    1326         xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45);
    1327 
    1328         /* Start  looping on the number of gaussian points: */
    1329         gauss=new GaussTria(2);
    1330         for(ig=gauss->begin();ig<gauss->end();ig++){
    1331 
    1332                 gauss->GaussPoint(ig);
    1333 
    1334                 tria->GetNodalFunctions(l1l2l3_tria,gauss);
    1335 
    1336                 //Get the coordinates of gauss point for each triangle in penta/quad
    1337                 r_tria=gauss->coord2-gauss->coord1;
    1338                 s_tria=-3/sqrt((double)3)*(gauss->coord1+gauss->coord2-2/3);
    1339 
    1340                 //Coordinates of gauss points in the reference triangle
    1341                 r_quad[0]=r_tria;
    1342                 s_quad[0]=1/sqrt((double)3)*s_tria-2/3;
    1343                 r_quad[1]=-1/sqrt((double)3)*s_tria+2/3;
    1344                 s_quad[1]=r_tria;
    1345                 r_quad[2]=-r_tria;
    1346                 s_quad[2]=-1/sqrt((double)3)*s_tria+2/3;
    1347                 r_quad[3]=1/sqrt((double)3)*s_tria-2/3;
    1348                 s_quad[3]=-r_tria;
    1349 
    1350                 //Get the nodal function of the quad for the gauss points of each triangle
    1351                 for(i=0;i<4;i++){
    1352                         l1l4_tria[i][0]=1.0/4.0*(
    1353                                         (r_quad[i]-1.0)*(s_quad[i]-1.0)
    1354                                         );
    1355                        
    1356                         l1l4_tria[i][1]=1.0/4.0*(
    1357                                          -(r_quad[i]+1.0)*(s_quad[i]-1.0)
    1358                                         );
    1359 
    1360                         l1l4_tria[i][2]=1.0/4.0*(
    1361                                         (r_quad[i]+1.0)*(s_quad[i]+1.0)
    1362                                         );
    1363                        
    1364                         l1l4_tria[i][3]=1.0/4.0*(
    1365                                          -(r_quad[i]-1.0)*(s_quad[i]+1.0)
    1366                                         );
    1367 
    1368                 } //for(i=0;i<4;i++)
    1369                
    1370                
    1371                 //Compute jacobian of each triangle
    1372                 for (i=0;i<4;i++){
    1373                         double complete_list[3][3]; //a third coordinate is needed for the jacobian determinant calculation, here it is zero
    1374                         for(j=0;j<3;j++){
    1375                                 complete_list[j][0]=xyz_tria[3*i+j][0];
    1376                                 complete_list[j][1]=xyz_tria[3*i+j][1];
    1377                                 complete_list[j][2]=0;
    1378                         }
    1379                         tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss);
    1380                 }
    1381 
    1382                 //Calculation of the z coordinate for the gaussian point ig for each triangle
    1383                 z_g[0]=z[0]*l1l2l3_tria[0]+z[1]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1384                 z_g[1]=z[1]*l1l2l3_tria[0]+z[2]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1385                 z_g[2]=z[2]*l1l2l3_tria[0]+z[3]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1386                 z_g[3]=z[3]*l1l2l3_tria[0]+z[0]*l1l2l3_tria[1]+z[4]*l1l2l3_tria[2];
    1387                
    1388                 //Loop on the triangles
    1389                 for(i=0;i<4;i++){
    1390 
    1391                         //Loop on the grids of the quad
    1392                         //Calculate the ice pressure
    1393                         air_pressure_tria=0;
    1394 
    1395                         //Now deal with water pressure:
    1396                         if(fill==WaterEnum){ //icefront ends in water
    1397                                 water_level_above_g_tria=min(0,z_g[i]);//0 if the gaussian point is above water level
    1398                                 water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
    1399                         }
    1400                         else if(fill==AirEnum){
    1401                                 water_pressure_tria=0;
    1402                         }
    1403                         else{
    1404                                 ISSMERROR("QuadPressureLoadStokes error message: unknow fill type for icefront boundary condition");
    1405                         }
    1406 
    1407                         //Add pressure
    1408                         pressure_tria = water_pressure_tria + air_pressure_tria;
    1409 
    1410                         pe_g[0]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nx[i];
    1411                         pe_g[1]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*ny[i];
    1412                         pe_g[2]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][0]*nz[i];
    1413                         pe_g[3]+= 0;
    1414                         pe_g[4]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nx[i];
    1415                         pe_g[5]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*ny[i];
    1416                         pe_g[6]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][1]*nz[i];
    1417                         pe_g[7]+= 0;
    1418                         pe_g[8]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nx[i];
    1419                         pe_g[9]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*ny[i];
    1420                         pe_g[10]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][2]*nz[i];
    1421                         pe_g[11]+= 0;
    1422                         pe_g[12]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nx[i];
    1423                         pe_g[13]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*ny[i];
    1424                         pe_g[14]+= J[i]*gauss->weight* pressure_tria*l1l4_tria[i][3]*nz[i];
    1425                         pe_g[15]+= 0;
    1426 
    1427                 } //for(i=0;i<4;i++)
    1428         }
    1429 
    1430         delete tria;
    1431         delete gauss;
    1432 }
    1433 /*}}}*/
    1434 /*FUNCTION Icefront::GetNormal {{{1*/
    1435 void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){
     874/*FUNCTION Icefront::GetSegmentNormal {{{1*/
     875void Icefront:: GetSegmentNormal(double* normal,double xyz_list[4][3]){
    1436876
    1437877        /*Build unit outward pointing vector*/
    … …  
    1447887        normal[0]= + vector[1]/norm;
    1448888        normal[1]= - vector[0]/norm;
     889}
     890/*}}}*/
     891/*FUNCTION Icefront::GetQuadNormal {{{1*/
     892void Icefront:: GetQuadNormal(double* normal,double xyz_list[4][3]){
     893
     894        /*Build unit outward pointing vector*/
     895        double AB[3];
     896        double AC[3];
     897        double norm;
     898
     899        AB[0]=xyz_list[1][0] - xyz_list[0][0];
     900        AB[1]=xyz_list[1][1] - xyz_list[0][1];
     901        AB[2]=xyz_list[1][2] - xyz_list[0][2];
     902        AC[0]=xyz_list[2][0] - xyz_list[0][0];
     903        AC[1]=xyz_list[2][1] - xyz_list[0][1];
     904        AC[2]=xyz_list[2][2] - xyz_list[0][2];
     905
     906        cross(normal,AB,AC);
     907        norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));
     908
     909        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
    1449910}
    1450911/*}}}*/

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

    r5772 r5808  
    8585                int*  GetExternalDofList(int approximation_enum,int setenum);
    8686                int   GetNumberOfDofs(int approximation_enum,int setenum);
    87                 void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    88                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
    89                 void  QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    90                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
    91                 void GetNormal(double* normal,double xyz_list[2][3]);
     87                void GetSegmentNormal(double* normal,double xyz_list[2][3]);
     88                void GetQuadNormal(double* normal,double xyz_list[4][3]);
    9289                ElementVector* NewElementVector(int approximation);
    9390                /*}}}*/

  • issm/trunk/src/c/shared/Numerics/GaussPoints.cpp

    r5744 r5808  
    17201720        GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
    17211721}/*}}}1*/
    1722 /*FUNCTION gaussPenta{{{1*/
    1723 void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus ) {
    1724         /*Gauss quadrature points for the pentahedron.*/
    1725 
    1726         /*  get the gauss points using the product of triangular and line rules  */
    1727         GaussLegendreTria(pngaus, pl1, pl2, pl3, pwgt, iord);
    1728         GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
    1729        
    1730 }/*}}}1*/

  • issm/trunk/src/c/shared/Numerics/GaussPoints.h

    r5744 r5808  
    1919void gaussQuad( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, int nigaus, int njgaus );
    2020void gaussHexa( double** pxgaus, double** pxwgt, double** pegaus, double** pewgt, double** pzgaus, double** pzwgt, int nigaus, int njgaus, int nkgaus );
    21 void gaussPenta( int* pngaus, double** pl1, double** pl2, double** pl3, double** pwgt, double** pzgaus, double** pzwgt, int iord, int nkgaus );
    2221
    2322#endif
Note: See TracChangeset for help on using the changeset viewer.