Changeset 20666

Timestamp:

05/30/16 12:00:16 (9 years ago)

Author:

Mathieu Morlighem

Message:

CHG: use diffusion equation to extrude from top

Location:

issm/trunk-jpl/src/c/analyses

Files:

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

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

@@ -49 +49 @@
 ElementMatrix* ExtrudeFromTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
-        /*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 */
-        int         dim;
         IssmDouble  Jdet,D;
         IssmDouble *xyz_list = NULL;
 
         /*Get dimension*/
+        int dim;
         element->FindParam(&dim,DomainDimensionEnum);
 
@@ -76 +62 @@
         /*Initialize Element vector and other vectors*/
         ElementMatrix* Ke     = element->NewElementMatrix();
-        IssmDouble*    B      = xNew<IssmDouble>(numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
 
         /*Retrieve all inputs and parameters*/
@@ -88 +73 @@
 
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                element->NodalFunctions(Bprime,gauss);
-                GetB(B,element,dim,xyz_list,gauss);
                 D=gauss->weight*Jdet;
+                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
-                TripleMultiply(B,1,numnodes,1,
-                                        &D,1,1,0,
-                                        Bprime,1,numnodes,0,
-                                        &Ke->values[0],1);
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
+                        }
+                }
         } 
 
         /*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 */
-        int         dim;
-        IssmDouble  Jdet,D,normal[2];
-        IssmDouble *xyz_list_top = NULL;
-
-        /*Get dimension*/
-        element->FindParam(&dim,DomainDimensionEnum);
-
-        /*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[dim-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->IsOnBase()) return NULL;
-
-        /*Intermediaries */
-        int         dim;
-        IssmDouble  Jdet,D,normal[3];
-        IssmDouble *xyz_list_base = NULL;
-
-        /*Get dimension*/
-        element->FindParam(&dim,DomainDimensionEnum);
-
-        /*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[dim-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){/*{{{*/
         return NULL;
 }/*}}}*/
-void           ExtrudeFromTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,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>(dim*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[i] = dbasis[(dim-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* CreateKMatrixBed(Element* element);
-                ElementMatrix* CreateKMatrixSurface(Element* element);
-                ElementMatrix* CreateKMatrixVolume(Element* element);
                 ElementVector* CreatePVector(Element* element);
-                void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
                 void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);