Changeset 16890


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

NEW: FS flowband

Location:
issm/trunk-jpl/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r16863 r16890  
    2727/*Finite Element Analysis*/
    2828ElementMatrix* AdjointBalancethicknessAnalysis::CreateKMatrix(Element* element){/*{{{*/
    29         _error_("not implemented yet");
     29
     30        BalancethicknessAnalysis* analysis = new BalancethicknessAnalysis();
     31        ElementMatrix* Ke = analysis->CreateKMatrix(element);
     32        delete analysis;
     33
     34        /*Transpose and return Ke*/
     35        Ke->Transpose();
     36        return Ke;
    3037}/*}}}*/
    3138ElementVector* AdjointBalancethicknessAnalysis::CreatePVector(Element* element){/*{{{*/

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

    r16862 r16890  
    3838/*Finite Element Analysis*/
    3939ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
    40         _error_("not implemented yet");
    41 }/*}}}*/
     40
     41        /*compute all stiffness matrices for this element*/
     42        ElementMatrix* Ke1=CreateKMatrixVolume(element);
     43        ElementMatrix* Ke2=CreateKMatrixSurface(element);
     44        ElementMatrix* Ke3=CreateKMatrixBed(element);
     45        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
     46
     47        /*clean-up and return*/
     48        delete Ke1;
     49        delete Ke2;
     50        delete Ke3;
     51        return Ke;
     52}/*}}}*/
     53ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
     54
     55        /*Intermediaries */
     56        IssmDouble  Jdet,D;
     57        IssmDouble *xyz_list = NULL;
     58
     59        /*Fetch number of nodes and dof for this finite element*/
     60        int numnodes = element->GetNumberOfNodes();
     61
     62        /*Initialize Element vector and other vectors*/
     63        ElementMatrix* Ke     = element->NewElementMatrix();
     64        IssmDouble*    B      = xNew<IssmDouble>(numnodes);
     65        IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
     66
     67        /*Retrieve all inputs and parameters*/
     68        element->GetVerticesCoordinates(&xyz_list);
     69
     70        /* Start  looping on the number of gaussian points: */
     71        Gauss* gauss=element->NewGauss(2);
     72        for(int ig=gauss->begin();ig<gauss->end();ig++){
     73                gauss->GaussPoint(ig);
     74
     75                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     76                element->NodalFunctions(Bprime,gauss);
     77                GetB(B,element,xyz_list,gauss);
     78                D=gauss->weight*Jdet;
     79
     80                TripleMultiply(B,1,numnodes,1,
     81                                        &D,1,1,0,
     82                                        Bprime,1,numnodes,0,
     83                                        &Ke->values[0],1);
     84        }
     85
     86        /*Clean up and return*/
     87        xDelete<IssmDouble>(xyz_list);
     88        xDelete<IssmDouble>(B);
     89        xDelete<IssmDouble>(Bprime);
     90        delete gauss;
     91        return Ke;
     92}
     93/*}}}*/
     94ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
     95
     96        if(!element->IsOnSurface()) return NULL;
     97
     98        /*Intermediaries */
     99        IssmDouble  Jdet,D,normal[2];
     100        IssmDouble *xyz_list_top = NULL;
     101
     102        /*Fetch number of nodes and dof for this finite element*/
     103        int numnodes = element->GetNumberOfNodes();
     104
     105        /*Initialize Element vector and other vectors*/
     106        ElementMatrix* Ke    = element->NewElementMatrix();
     107        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     108
     109        /*Retrieve all inputs and parameters*/
     110        element->GetVerticesCoordinatesTop(&xyz_list_top);
     111        element->NormalTop(&normal[0],xyz_list_top);
     112
     113        /* Start  looping on the number of gaussian points: */
     114        Gauss* gauss=element->NewGaussTop(2);
     115        for(int ig=gauss->begin();ig<gauss->end();ig++){
     116                gauss->GaussPoint(ig);
     117
     118                element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
     119                element->NodalFunctions(basis,gauss);
     120                D = - gauss->weight*Jdet*normal[1];
     121
     122                TripleMultiply(basis,1,numnodes,1,
     123                                        &D,1,1,0,
     124                                        basis,1,numnodes,0,
     125                                        &Ke->values[0],1);
     126        }
     127
     128        /*Clean up and return*/
     129        xDelete<IssmDouble>(xyz_list_top);
     130        xDelete<IssmDouble>(basis);
     131        delete gauss;
     132        return Ke;
     133}
     134/*}}}*/
     135ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
     136
     137        if(!element->IsOnBed()) return NULL;
     138
     139        /*Intermediaries */
     140        IssmDouble  Jdet,D,normal[2];
     141        IssmDouble *xyz_list_base = NULL;
     142
     143        /*Fetch number of nodes and dof for this finite element*/
     144        int numnodes = element->GetNumberOfNodes();
     145
     146        /*Initialize Element vector and other vectors*/
     147        ElementMatrix* Ke    = element->NewElementMatrix();
     148        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     149
     150        /*Retrieve all inputs and parameters*/
     151        element->GetVerticesCoordinatesBase(&xyz_list_base);
     152        element->NormalBase(&normal[0],xyz_list_base);
     153
     154        /* Start  looping on the number of gaussian points: */
     155        Gauss* gauss=element->NewGaussBase(2);
     156        for(int ig=gauss->begin();ig<gauss->end();ig++){
     157                gauss->GaussPoint(ig);
     158
     159                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     160                element->NodalFunctions(basis,gauss);
     161                D = - gauss->weight*Jdet*normal[1];
     162
     163                TripleMultiply(basis,1,numnodes,1,
     164                                        &D,1,1,0,
     165                                        basis,1,numnodes,0,
     166                                        &Ke->values[0],1);
     167        }
     168
     169        /*Clean up and return*/
     170        xDelete<IssmDouble>(xyz_list_base);
     171        xDelete<IssmDouble>(basis);
     172        delete gauss;
     173        return Ke;
     174}
     175/*}}}*/
    42176ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/
    43177        return NULL;
    44178}/*}}}*/
     179void ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     180        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
     181                where hi is the interpolation function for node i.*/
     182
     183        /*Fetch number of nodes for this finite element*/
     184        int numnodes = element->GetNumberOfNodes();
     185
     186        /*Get nodal functions derivatives*/
     187        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     188        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     189
     190        /*Build B: */
     191        for(int i=0;i<numnodes;i++){
     192                B[i] = dbasis[1*numnodes+i]; 
     193        }
     194
     195        /*Clean-up*/
     196        xDelete<IssmDouble>(dbasis);
     197}
     198/*}}}*/
    45199void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    46200           _error_("not implemented yet");

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h

    r16782 r16890  
    2222                /*Finite element Analysis*/
    2323                ElementMatrix* CreateKMatrix(Element* element);
     24                ElementMatrix* CreateKMatrixVolume(Element* element);
     25                ElementMatrix* CreateKMatrixSurface(Element* element);
     26                ElementMatrix* CreateKMatrixBed(Element* element);
    2427                ElementVector* CreatePVector(Element* element);
     28                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    2529                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2630                void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

    r16862 r16890  
    6666                case Mesh2DhorizontalEnum:
    6767                        basalelement = element;
     68                        break;
     69                case Mesh2DverticalEnum:
     70                        if(!element->IsOnBed()) return NULL;
     71                        basalelement = element->SpawnBasalElement();
    6872                        break;
    6973                case Mesh3DEnum:

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16876 r16890  
    12131213
    12141214}/*}}}*/
     1215/*FUNCTION Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){{{*/
     1216void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){
     1217
     1218        int        indices[2];
     1219        IssmDouble xyz_list[NUMVERTICES][3];
     1220
     1221        /*Element XYZ list*/
     1222        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     1223
     1224        /*Allocate Output*/
     1225        IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
     1226        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
     1227        for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
     1228
     1229        /*Assign output pointer*/
     1230        *pxyz_list = xyz_list_edge;
     1231
     1232}/*}}}*/
    12151233/*FUNCTION Tria::NormalSection{{{*/
    12161234void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
    … …  
    20862104}
    20872105/*}}}*/
     2106/*FUNCTION Tria::JacobianDeterminantTop{{{*/
     2107void Tria::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){
     2108
     2109        _assert_(gauss->Enum()==GaussTriaEnum);
     2110        this->GetSegmentJacobianDeterminant(pJdet,xyz_list_top,(GaussTria*)gauss);
     2111
     2112}
     2113/*}}}*/
    20882114/*FUNCTION Tria::HasEdgeOnBed {{{*/
    20892115bool Tria::HasEdgeOnBed(){
    … …  
    23102336        int indices[2];
    23112337        this->EdgeOnBedIndices(&indices[0],&indices[1]);
     2338        return new GaussTria(indices[0],indices[1],order);
     2339}
     2340/*}}}*/
     2341/*FUNCTION Tria::NewGaussTop(int order){{{*/
     2342Gauss* Tria::NewGaussTop(int order){
     2343
     2344        int indices[2];
     2345        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
    23122346        return new GaussTria(indices[0],indices[1],order);
    23132347}
    … …  
    24122446/*FUNCTION Tria::NormalBase {{{*/
    24132447void Tria::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
     2448
     2449        /*Build unit outward pointing vector*/
     2450        IssmDouble vector[2];
     2451        IssmDouble norm;
     2452
     2453        vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
     2454        vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
     2455
     2456        norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
     2457
     2458        bed_normal[0]= + vector[1]/norm;
     2459        bed_normal[1]= - vector[0]/norm;
     2460}
     2461/*}}}*/
     2462/*FUNCTION Tria::NormalTop {{{*/
     2463void Tria::NormalTop(IssmDouble* bed_normal,IssmDouble* xyz_list){
    24142464
    24152465        /*Build unit outward pointing vector*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16886 r16890  
    121121                void        GetVerticesCoordinates(IssmDouble** pxyz_list);
    122122                void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
    123                 void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
     123                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
    124124                void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
    125125                void        InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum);
    … …  
    270270                IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
    271271                void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
    272                 void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
     272                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
    273273                void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
    274274                IssmDouble     GetMaterialParameter(int enum_in);
    … …  
    294294                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    295295                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    296                 void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
     296                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    297297                IssmDouble     MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
    298298                Gauss*         NewGauss(void);
    … …  
    303303                Gauss*         NewGaussBase(int order);
    304304                Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
    305                 Gauss*         NewGaussTop(int order){_error_("not implemented yet");};
     305                Gauss*         NewGaussTop(int order);
    306306                ElementVector* NewElementVector(int approximation_enum);
    307307                ElementMatrix* NewElementMatrix(int approximation_enum);
Note: See TracChangeset for help on using the changeset viewer.