Ignore:
Timestamp:
11/01/19 12:01:57 (5 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 24310

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        55*.obj
        66*.exe
         7*.mod
        78appscan.*
        89ouncemake*
        … …  
        1415probe.results
        1516stXXXX*
        16 *.deps
        17 *.dirstamp
         17.deps
         18.dirstamp
        1819.libs
        1920issm
        … …  
        2122issm_slr
        2223issm_ocean
        23 lnb_param.mod
        24 lovenb_sub.mod
        25 model.mod
        26 util.mod
         24issm_dakota

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

    r23189 r24313  
    3333/*Reference Element numerics*/
    3434void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    35         /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
     35        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
    3636         * point specified by gauss_basis:
    3737         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    … …  
    8989/*}}}*/
    9090void TriaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    91         /*The Jacobian is constant over the element, discard the gaussian points. 
     91        /*The Jacobian is constant over the element, discard the gaussian points.
    9292         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    9393
    … …  
    106106/*}}}*/
    107107void TriaRef::GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    108         /*The Jacobian determinant is constant over the element, discard the gaussian points. 
     108        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    109109         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    110110        IssmDouble J[2][2];
    … …  
    120120/*}}}*/
    121121void TriaRef::GetJacobianDeterminant3D(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    122         /*The Jacobian determinant is constant over the element, discard the gaussian points. 
     122        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    123123         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    124124        IssmDouble J[2][2];
    … …  
    157157                case NoneEnum:
    158158                        return;
    159                 case P0Enum:
     159                case P0Enum: case P0DGEnum:
    160160                        basis[0]=1.;
    161161                        return;
    … …  
    202202void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, Gauss* gauss,int finiteelement){/*{{{*/
    203203
    204         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
     204        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    205205         * actual coordinate system): */
    206206        IssmDouble    Jinv[2][2];
    … …  
    211211        /*Get nodal functions derivatives in reference triangle*/
    212212        IssmDouble dbasis_ref[2*NUMNODESMAX];
    213         GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement); 
     213        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss,finiteelement);
    214214
    215215        /*Get Jacobian invert: */
    216216        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    217217
    218         /*Build dbasis: 
     218        /*Build dbasis:
    219219         * [dhi/dx]= Jinv*[dhi/dr]
    220220         * [dhi/dy]       [dhi/ds]
    … …  
    227227}
    228228/*}}}*/
    229 void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    230         /*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
    231          *
    232          * and phi1=phi3 phi2=phi4
    233          *
    234          * We assume B has been allocated already, of size: 1x4
    235          */
    236 
    237         /*Fetch number of nodes for this finite element*/
    238         int numnodes = this->NumberofNodes(finiteelement);
    239 
    240         /*Get nodal functions*/
    241         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    242         GetNodalFunctions(basis,gauss,finiteelement);
    243 
    244         /*Build B for this segment*/
    245         B[0] = +basis[index1];
    246         B[1] = +basis[index2];
    247         B[2] = -basis[index1];
    248         B[3] = -basis[index2];
    249 
    250         /*Clean-up*/
    251         xDelete<IssmDouble>(basis);
    252 }
    253 /*}}}*/
    254 void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,Gauss* gauss, int index1,int index2,int finiteelement){/*{{{*/
    255         /*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
    256          *
    257          * and phi1=phi3 phi2=phi4
    258          *
    259          * We assume Bprime has been allocated already, of size: 1x4
    260          */
    261 
    262         /*Fetch number of nodes for this finite element*/
    263         int numnodes = this->NumberofNodes(finiteelement);
    264 
    265         /*Get nodal functions*/
    266         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    267         GetNodalFunctions(basis,gauss,finiteelement);
    268 
    269         /*Build B'*/
    270         Bprime[0] = basis[index1];
    271         Bprime[1] = basis[index2];
    272         Bprime[2] = basis[index1];
    273         Bprime[3] = basis[index2];
    274 
    275         /*Clean-up*/
    276         xDelete<IssmDouble>(basis);
    277 }
    278 /*}}}*/
    279229void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    280         /*The Jacobian determinant is constant over the element, discard the gaussian points. 
     230        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    281231         * J is assumed to have been allocated*/
    282232
    … …  
    305255
    306256        switch(finiteelement){
     257                case P0Enum: case P0DGEnum:
     258                        basis[0]=triabasis[0];
     259                        xDelete<IssmDouble>(triabasis);
     260                        return;
    307261                case P1Enum: case P1DGEnum:
    308262                        basis[0]=triabasis[index1];
    … …  
    330284}
    331285/*}}}*/
     286void TriaRef::GetSegmentNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list_tria,Gauss* gauss,int index1,int index2,int finiteelement){/*{{{*/
     287        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     288
     289        _assert_(index1>=0 && index1<3);
     290        _assert_(index2>=0 && index2<3);
     291
     292        /*Fetch number of nodes for this finite element*/
     293        int numnodes = this->NumberofNodes(finiteelement);
     294
     295        /*Get nodal functions*/
     296        IssmDouble* dtriabasis=xNew<IssmDouble>(2*numnodes);
     297        GetNodalFunctionsDerivatives(dtriabasis,xyz_list_tria,gauss,finiteelement);
     298
     299        switch(finiteelement){
     300                case P1Enum: case P1DGEnum:
     301                        dbasis[2*0+0] = dtriabasis[numnodes*0+index1];
     302                        dbasis[2*0+1] = dtriabasis[numnodes*1+index1];
     303                        dbasis[2*1+0] = dtriabasis[numnodes*0+index2];
     304                        dbasis[2*1+1] = dtriabasis[numnodes*1+index2];
     305                        break;
     306                default:
     307                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
     308        }
     309
     310        /*Clean up*/
     311        xDelete<IssmDouble>(dtriabasis);
     312}
     313/*}}}*/
    332314void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,Gauss* gauss_in,int finiteelement){/*{{{*/
    333         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
     315        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    334316         * natural coordinate system) at the gaussian point. */
    335317
    … …  
    341323
    342324        switch(finiteelement){
    343                 case P0Enum:
     325                case P0Enum: case P0DGEnum:
    344326                        /*Nodal function 1*/
    345327                        dbasis[NUMNODESP0*0+0] = 0.;
    … …  
    484466                case NoneEnum:                return 0;
    485467                case P0Enum:                  return NUMNODESP0;
     468                case P0DGEnum:                return NUMNODESP0;
    486469                case P1Enum:                  return NUMNODESP1;
    487470                case P1DGEnum:                return NUMNODESP1;
Note: See TracChangeset for help on using the changeset viewer.