Ignore:
Timestamp:
07/22/10 15:46:39 (15 years ago)
Author:
Mathieu Morlighem
Message:

moved all common routines of Tria and TriaVertexInput to TriaRef

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r4546 r4769  
    2424/*}}}*/
    2525/*FUNCTION TriaVertexInput::TriaVertexInput(int in_enum_type,double* values){{{1*/
    26 TriaVertexInput::TriaVertexInput(int in_enum_type,double* in_values){
    27 
     26TriaVertexInput::TriaVertexInput(int in_enum_type,double* in_values)
     27        :TriaRef(1)
     28{
     29
     30        /*Set TriaRef*/
     31        this->SetElementType(P1Enum,0);
     32        this->element_type=P1Enum;
     33
     34        /*Set Enum*/
    2835        enum_type=in_enum_type;
     36
     37        /*Set values*/
    2938        values[0]=in_values[0];
    3039        values[1]=in_values[1];
     
    311320
    312321/*Intermediary*/
    313 /*FUNCTION TriaVertexInput::GetNodalFunctions {{{1*/
    314 void TriaVertexInput::GetNodalFunctions(double* l1l2l3, double* gauss){
    315        
    316         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    317 
    318         /*First nodal function: */
    319         l1l2l3[0]=gauss[0];
    320 
    321         /*Second nodal function: */
    322         l1l2l3[1]=gauss[1];
    323 
    324         /*Third nodal function: */
    325         l1l2l3[2]=gauss[2];
    326 
    327 }
    328 /*}}}*/
    329 /*FUNCTION TriaVertexInput::GetB {{{1*/
    330 void TriaVertexInput::GetB(double* B, double* xyz_list, double* gauss){
    331         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    332          * For grid i, Bi can be expressed in the actual coordinate system
    333          * by:
    334          *       Bi=[ dh/dx           0    ]
    335          *          [   0           dh/dy  ]
    336          *          [ 1/2*dh/dy  1/2*dh/dx ]
    337          * where h is the interpolation function for grid i.
    338          *
    339          * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
    340          */
    341 
    342         int i;
    343         const int NDOF2=2;
    344         const int numgrids=3;
    345 
    346         double dh1dh3[NDOF2][numgrids];
    347 
    348 
    349         /*Get dh1dh2dh3 in actual coordinate system: */
    350         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
    351 
    352         /*Build B: */
    353         for (i=0;i<numgrids;i++){
    354                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
    355                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
    356                 *(B+NDOF2*numgrids*1+NDOF2*i)=0;
    357                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
    358                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];
    359                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];
    360         }
    361 }
    362 /*}}}*/
    363 /*FUNCTION TriaVertexInput::GetNodalFunctionsDerivatives {{{1*/
    364 void TriaVertexInput::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){
    365 
    366         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    367          * actual coordinate system: */
    368         int i;
    369         const int NDOF2=2;
    370         const int numgrids=3;
    371         double dh1dh3_ref[NDOF2][numgrids];
    372         double Jinv[NDOF2][NDOF2];
    373 
    374         /*Get derivative values with respect to parametric coordinate system: */
    375         GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss);
    376 
    377         /*Get Jacobian invert: */
    378         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    379 
    380         /*Build dh1dh3:
    381          *
    382          * [dhi/dx]= Jinv*[dhi/dr]
    383          * [dhi/dy]       [dhi/ds]
    384          */
    385         for (i=0;i<numgrids;i++){
    386                 dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
    387                 dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
    388         }
    389 
    390 }
    391 /*}}}*/
    392 /*FUNCTION TriaVertexInput::GetNodalFunctionsDerivativesReference {{{1*/
    393 void TriaVertexInput::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){
    394 
    395         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    396          * natural coordinate system) at the gaussian point. */
    397 
    398         const int NDOF2=2;
    399         const int numgrids=3;
    400 
    401         /*First nodal function: */
    402         *(dl1dl3+numgrids*0+0)=-0.5;
    403         *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
    404 
    405         /*Second nodal function: */
    406         *(dl1dl3+numgrids*0+1)=0.5;
    407         *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
    408 
    409         /*Third nodal function: */
    410         *(dl1dl3+numgrids*0+2)=0;
    411         *(dl1dl3+numgrids*1+2)=1.0/SQRT3;
    412 
    413 }
    414 /*}}}*/
    415 /*FUNCTION TriaVertexInput::GetJacobian {{{1*/
    416 void TriaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss){
    417 
    418         /*The Jacobian is constant over the element, discard the gaussian points.
    419          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    420 
    421         const int NDOF2=2;
    422         const int numgrids=3;
    423         double x1,y1,x2,y2,x3,y3;
    424 
    425         x1=*(xyz_list+numgrids*0+0);
    426         y1=*(xyz_list+numgrids*0+1);
    427         x2=*(xyz_list+numgrids*1+0);
    428         y2=*(xyz_list+numgrids*1+1);
    429         x3=*(xyz_list+numgrids*2+0);
    430         y3=*(xyz_list+numgrids*2+1);
    431 
    432 
    433         *(J+NDOF2*0+0)=0.5*(x2-x1);
    434         *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
    435         *(J+NDOF2*0+1)=0.5*(y2-y1);
    436         *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
    437 }
    438 /*}}}*/
    439 /*FUNCTION TriaVertexInput::GetJacobianInvert {{{1*/
    440 void TriaVertexInput::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
    441 
    442         double Jdet;
    443         const int NDOF2=2;
    444         const int numgrids=3;
    445 
    446         /*Call Jacobian routine to get the jacobian:*/
    447         GetJacobian(Jinv, xyz_list, gauss);
    448 
    449         /*Invert Jacobian matrix: */
    450         MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
    451 
    452 }
    453 /*}}}*/
    454322/*FUNCTION TriaVertexInput::SquareMin{{{1*/
    455323void TriaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
Note: See TracChangeset for help on using the changeset viewer.