Changeset 16890

Timestamp:

11/22/13 11:10:22 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: FS flowband

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointBalancethicknessAnalysis.cpp (modified) (1 diff)
• analyses/ExtrudeFromBaseAnalysis.cpp (modified) (1 diff)
• analyses/ExtrudeFromBaseAnalysis.h (modified) (1 diff)
• analyses/L2ProjectionBaseAnalysis.cpp (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (4 diffs)
• classes/Elements/Tria.h (modified) (4 diffs)

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

@@ -27 +27 @@
 /*Finite Element Analysis*/
 ElementMatrix* AdjointBalancethicknessAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
+
+        BalancethicknessAnalysis* analysis = new BalancethicknessAnalysis();
+        ElementMatrix* Ke = analysis->CreateKMatrix(element);
+        delete analysis;
+
+        /*Transpose and return Ke*/
+        Ke->Transpose();
+        return Ke;
 }/*}}}*/
 ElementVector* AdjointBalancethicknessAnalysis::CreatePVector(Element* element){/*{{{*/

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

@@ -38 +38 @@
 /*Finite Element Analysis*/
 ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
-        _error_("not implemented yet");
-}/*}}}*/
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=CreateKMatrixVolume(element);
+        ElementMatrix* Ke2=CreateKMatrixSurface(element);
+        ElementMatrix* Ke3=CreateKMatrixBed(element);
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
+
+        /*clean-up and return*/
+        delete Ke1;
+        delete Ke2;
+        delete Ke3;
+        return Ke;
+}/*}}}*/
+ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble  Jdet,D;
+        IssmDouble *xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(Bprime,gauss);
+                GetB(B,element,xyz_list,gauss);
+                D=gauss->weight*Jdet;
+
+                TripleMultiply(B,1,numnodes,1,
+                                        &D,1,1,0,
+                                        Bprime,1,numnodes,0,
+                                        &Ke->values[0],1);
+        } 
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
+ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
+
+        if(!element->IsOnSurface()) return NULL;
+
+        /*Intermediaries */
+        IssmDouble  Jdet,D,normal[2];
+        IssmDouble *xyz_list_top = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke    = element->NewElementMatrix();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinatesTop(&xyz_list_top);
+        element->NormalTop(&normal[0],xyz_list_top);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGaussTop(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
+                element->NodalFunctions(basis,gauss);
+                D = - gauss->weight*Jdet*normal[1]; 
+
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Ke->values[0],1);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list_top);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
+ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
+
+        if(!element->IsOnBed()) return NULL;
+
+        /*Intermediaries */
+        IssmDouble  Jdet,D,normal[2];
+        IssmDouble *xyz_list_base = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke    = element->NewElementMatrix();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinatesBase(&xyz_list_base);
+        element->NormalBase(&normal[0],xyz_list_base);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGaussBase(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+                element->NodalFunctions(basis,gauss);
+                D = - gauss->weight*Jdet*normal[1]; 
+
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Ke->values[0],1);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list_base);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        return Ke;
+}
+/*}}}*/
 ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/
         return NULL;
 }/*}}}*/
+void ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
+                where hi is the interpolation function for node i.*/
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B[i] = dbasis[1*numnodes+i];  
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+}
+/*}}}*/
 void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
            _error_("not implemented yet");

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

@@ -22 +22 @@
                 /*Finite element Analysis*/
                 ElementMatrix* CreateKMatrix(Element* element);
+                ElementMatrix* CreateKMatrixVolume(Element* element);
+                ElementMatrix* CreateKMatrixSurface(Element* element);
+                ElementMatrix* CreateKMatrixBed(Element* element);
                 ElementVector* CreatePVector(Element* element);
+                void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
                 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);

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

@@ -66 +66 @@
                 case Mesh2DhorizontalEnum:
                         basalelement = element;
+                        break;
+                case Mesh2DverticalEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
                         break;
                 case Mesh3DEnum:

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

@@ -1213 +1213 @@
 
 }/*}}}*/
+/*FUNCTION Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){{{*/
+void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){
+
+        int        indices[2];
+        IssmDouble xyz_list[NUMVERTICES][3];
+
+        /*Element XYZ list*/
+        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+
+        /*Allocate Output*/
+        IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
+        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
+        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];
+
+        /*Assign output pointer*/
+        *pxyz_list = xyz_list_edge;
+
+}/*}}}*/
 /*FUNCTION Tria::NormalSection{{{*/
 void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){
@@ -2086 +2104 @@
 }
 /*}}}*/
+/*FUNCTION Tria::JacobianDeterminantTop{{{*/
+void Tria::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){
+
+        _assert_(gauss->Enum()==GaussTriaEnum);
+        this->GetSegmentJacobianDeterminant(pJdet,xyz_list_top,(GaussTria*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tria::HasEdgeOnBed {{{*/
 bool Tria::HasEdgeOnBed(){
@@ -2310 +2336 @@
         int indices[2];
         this->EdgeOnBedIndices(&indices[0],&indices[1]);
+        return new GaussTria(indices[0],indices[1],order);
+}
+/*}}}*/
+/*FUNCTION Tria::NewGaussTop(int order){{{*/
+Gauss* Tria::NewGaussTop(int order){
+
+        int indices[2];
+        this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
         return new GaussTria(indices[0],indices[1],order);
 }
@@ -2412 +2446 @@
 /*FUNCTION Tria::NormalBase {{{*/
 void Tria::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){
+
+        /*Build unit outward pointing vector*/
+        IssmDouble vector[2];
+        IssmDouble norm;
+
+        vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+        vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+
+        norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
+
+        bed_normal[0]= + vector[1]/norm;
+        bed_normal[1]= - vector[0]/norm;
+}
+/*}}}*/
+/*FUNCTION Tria::NormalTop {{{*/
+void Tria::NormalTop(IssmDouble* bed_normal,IssmDouble* xyz_list){
 
         /*Build unit outward pointing vector*/

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

@@ -121 +121 @@
                 void        GetVerticesCoordinates(IssmDouble** pxyz_list);
                 void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
-                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");};
+                void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
                 void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
                 void        InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum);
@@ -270 +270 @@
                 IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
                 void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
-                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
+                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
                 void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
                 IssmDouble     GetMaterialParameter(int enum_in);
@@ -294 +294 @@
                 void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
                 void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
-                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
+                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
                 IssmDouble     MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
                 Gauss*         NewGauss(void);
@@ -303 +303 @@
                 Gauss*         NewGaussBase(int order);
                 Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
-                Gauss*         NewGaussTop(int order){_error_("not implemented yet");};
+                Gauss*         NewGaussTop(int order);
                 ElementVector* NewElementVector(int approximation_enum);
                 ElementMatrix* NewElementMatrix(int approximation_enum);