Changeset 16889


Ignore:
Timestamp:
11/22/13 10:36:44 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: K matrix bal vel

Location:
issm/trunk-jpl/src/c/analyses
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp

    r16860 r16889  
    6060/*Finite Element Analysis*/
    6161ElementMatrix* BalancevelocityAnalysis::CreateKMatrix(Element* element){/*{{{*/
    62         _error_("not implemented yet");
     62
     63        /*Intermediaries */
     64        IssmDouble  dhdt,mb,ms,Jdet;
     65        IssmDouble  h,gamma,thickness;
     66        IssmDouble  hnx,hny,dhnx[2],dhny[2];
     67        IssmDouble *xyz_list = NULL;
     68
     69        /*Fetch number of nodes and dof for this finite element*/
     70        int numnodes = element->GetNumberOfNodes();
     71
     72        /*Initialize Element matrix and vectors*/
     73        ElementMatrix* Ke     = element->NewElementMatrix();
     74        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     75        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     76        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     77        IssmDouble*    HNx    = xNew<IssmDouble>(numnodes);
     78        IssmDouble*    HNy    = xNew<IssmDouble>(numnodes);
     79        IssmDouble*    H      = xNew<IssmDouble>(numnodes);
     80        IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
     81        IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
     82
     83        /*Retrieve all Inputs and parameters: */
     84        element->GetVerticesCoordinates(&xyz_list);
     85        Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input);
     86        h = element->CharacteristicLength();
     87
     88        /*Get vector N for all nodes and build HNx and HNy*/
     89        element->GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
     90        element->GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
     91        element->GetInputListOnNodes(H,ThicknessEnum);
     92        for(int i=0;i<numnodes;i++){
     93                IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
     94                HNx[i] = -H[i]*Nx[i]/norm;
     95                HNy[i] = -H[i]*Ny[i]/norm;
     96        }
     97
     98        /*Start looping on the number of gaussian points:*/
     99        Gauss* gauss=element->NewGauss(2);
     100        for(int ig=gauss->begin();ig<gauss->end();ig++){
     101                gauss->GaussPoint(ig);
     102
     103                H_input->GetInputValue(&thickness,gauss);
     104                if(thickness<50.) thickness=50.;
     105                element->ValueP1DerivativesOnGauss(&dhnx[0],HNx,xyz_list,gauss);
     106                element->ValueP1DerivativesOnGauss(&dhny[0],HNy,xyz_list,gauss);
     107                element->ValueP1OnGauss(&hnx,HNx,gauss);
     108                element->ValueP1OnGauss(&hny,HNy,gauss);
     109
     110                gamma=h/(2.*thickness+1.e-10);
     111
     112                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     113                element->NodalFunctions(basis,gauss);
     114                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     115
     116                for(int i=0;i<numnodes;i++){
     117                        for(int j=0;j<numnodes;j++){
     118                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
     119                                                        (basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*
     120                                                        (basis[j]*(dhnx[0]+dhny[1])  + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
     121                                                        );
     122                        }
     123                }
     124        }
     125
     126        /*Clean up and return*/
     127        xDelete<IssmDouble>(xyz_list);
     128        xDelete<IssmDouble>(basis);
     129        xDelete<IssmDouble>(dbasis);
     130        xDelete<IssmDouble>(H);
     131        xDelete<IssmDouble>(Nx);
     132        xDelete<IssmDouble>(Ny);
     133        xDelete<IssmDouble>(HNx);
     134        xDelete<IssmDouble>(HNy);
     135        xDelete<IssmDouble>(B);
     136        delete gauss;
     137        return Ke;
    63138}/*}}}*/
    64139ElementVector* BalancevelocityAnalysis::CreatePVector(Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeXAnalysis.cpp

    r16860 r16889  
    4646/*Finite Element Analysis*/
    4747ElementMatrix* SmoothedSurfaceSlopeXAnalysis::CreateKMatrix(Element* element){/*{{{*/
    48         _error_("not implemented yet");
     48
     49        /* Intermediaries */
     50        IssmDouble  Jdet,thickness,l=8.;
     51        IssmDouble *xyz_list = NULL;
     52
     53        /*Fetch number of nodes and dof for this finite element*/
     54        int numnodes = element->GetNumberOfNodes();
     55
     56        /*Initialize Element matrix and vectors*/
     57        ElementMatrix* Ke     = element->NewElementMatrix();
     58        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     59        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     60
     61        /*Retrieve all inputs and parameters*/
     62        element->GetVerticesCoordinates(&xyz_list);
     63        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     64
     65        /* Start looping on the number of gaussian points: */
     66        Gauss* gauss=element->NewGauss(2);
     67        for(int ig=gauss->begin();ig<gauss->end();ig++){
     68                gauss->GaussPoint(ig);
     69
     70                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     71                thickness_input->GetInputValue(&thickness,gauss);
     72                if(thickness<50.) thickness=50.;
     73
     74                element->NodalFunctions(basis,gauss);
     75                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     76
     77                for(int i=0;i<numnodes;i++){
     78                        for(int j=0;j<numnodes;j++){
     79                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
     80                                                        basis[i]*basis[j]
     81                                                        +(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
     82                                                        );
     83                        }
     84                }
     85        }
     86
     87        /*Clean up and return*/
     88        delete gauss;
     89        xDelete<IssmDouble>(dbasis);
     90        xDelete<IssmDouble>(basis);
     91        xDelete<IssmDouble>(xyz_list);
     92        return Ke;
    4993}/*}}}*/
    5094ElementVector* SmoothedSurfaceSlopeXAnalysis::CreatePVector(Element* element){/*{{{*/

  • issm/trunk-jpl/src/c/analyses/SmoothedSurfaceSlopeYAnalysis.cpp

    r16860 r16889  
    4646/*Finite Element Analysis*/
    4747ElementMatrix* SmoothedSurfaceSlopeYAnalysis::CreateKMatrix(Element* element){/*{{{*/
    48         _error_("not implemented yet");
     48
     49        /* Intermediaries */
     50        IssmDouble  Jdet,thickness,l=8.;
     51        IssmDouble *xyz_list = NULL;
     52
     53        /*Fetch number of nodes and dof for this finite element*/
     54        int numnodes = element->GetNumberOfNodes();
     55
     56        /*Initialize Element matrix and vectors*/
     57        ElementMatrix* Ke     = element->NewElementMatrix();
     58        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
     59        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     60
     61        /*Retrieve all inputs and parameters*/
     62        element->GetVerticesCoordinates(&xyz_list);
     63        Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     64
     65        /* Start looping on the number of gaussian points: */
     66        Gauss* gauss=element->NewGauss(2);
     67        for(int ig=gauss->begin();ig<gauss->end();ig++){
     68                gauss->GaussPoint(ig);
     69
     70                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     71                thickness_input->GetInputValue(&thickness,gauss);
     72                if(thickness<50.) thickness=50.;
     73
     74                element->NodalFunctions(basis,gauss);
     75                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     76
     77                for(int i=0;i<numnodes;i++){
     78                        for(int j=0;j<numnodes;j++){
     79                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
     80                                                        basis[i]*basis[j]
     81                                                        +(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
     82                                                        );
     83                        }
     84                }
     85        }
     86
     87        /*Clean up and return*/
     88        delete gauss;
     89        xDelete<IssmDouble>(dbasis);
     90        xDelete<IssmDouble>(basis);
     91        xDelete<IssmDouble>(xyz_list);
     92        return Ke;
    4993}/*}}}*/
    5094ElementVector* SmoothedSurfaceSlopeYAnalysis::CreatePVector(Element* element){/*{{{*/
Note: See TracChangeset for help on using the changeset viewer.