Changeset 17179

Timestamp:

01/28/14 08:06:20 (11 years ago)

Author:

Mathieu Morlighem

Message:

BUG: some fixes

Location:

issm/trunk-jpl/src/c/analyses

Files:

• ExtrudeFromTopAnalysis.cpp (modified) (1 diff)
• ExtrudeFromTopAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp

@@ -48 +48 @@
 }/*}}}*/
 ElementMatrix* ExtrudeFromTopAnalysis::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* ExtrudeFromTopAnalysis::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* ExtrudeFromTopAnalysis::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* ExtrudeFromTopAnalysis::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* ExtrudeFromTopAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
-}/*}}}*/
+        return NULL;
+}/*}}}*/
+void ExtrudeFromTopAnalysis::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 ExtrudeFromTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
            _error_("not implemented yet");

issm/trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h

@@ -25 +25 @@
                 ElementMatrix* CreateJacobianMatrix(Element* element);
                 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);