Changeset 16993 for issm/trunk-jpl/src/c/classes/Elements/Seg.cpp

Timestamp:

12/03/13 10:36:26 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: moving all Stiffness matrices and Jacobian matrices to analyses

File:

• issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Seg.cpp

@@ -138 +138 @@
 }
 /*}}}*/
-/*FUNCTION Seg::CreateKMatrixFreeSurfaceTop {{{*/
-ElementMatrix* Seg::CreateKMatrixFreeSurfaceTop(void){
-
-        /*Intermediaries */
-        int        stabilization;
-        IssmDouble Jdet,D_scalar,dt,h;
-        IssmDouble vx,vel;
-        IssmDouble xyz_list[NUMVERTICES][3];
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = this->NumberofNodes();
-
-        /*Initialize Element matrix and vectors*/
-        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
-        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
-
-        /*Retrieve all inputs and parameters*/
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
-        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-        h=this->GetSize();
-
-        /* Start  looping on the number of gaussian points: */
-        GaussSeg *gauss=new GaussSeg(2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-                GetNodalFunctions(basis,gauss);
-
-                vx_input->GetInputValue(&vx,gauss);
-
-                D_scalar=gauss->weight*Jdet;
-
-                TripleMultiply(basis,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        basis,1,numnodes,0,
-                                        &Ke->values[0],1);
-
-                GetNodalFunctions(B,gauss);
-                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
-
-                D_scalar=dt*gauss->weight*Jdet*vx;
-                TripleMultiply(B,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        Bprime,1,numnodes,0,
-                                        &Ke->values[0],1);
-
-                if(stabilization==2){
-                        /*Streamline upwinding*/
-                        vel=fabs(vx)+1.e-8;
-                        D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
-                }
-                else if(stabilization==1){
-                        /*SSA*/
-                        vx_input->GetInputAverage(&vx);
-                        D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
-                }
-                if(stabilization==1 || stabilization==2){
-                        TripleMultiply(Bprime,1,numnodes,1,
-                                                &D_scalar,1,1,0,
-                                                Bprime,1,numnodes,0,
-                                                &Ke->values[0],1);
-                }
-        }
-
-        /*Clean up and return*/
-        xDelete<IssmDouble>(basis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
-        delete gauss;
-        return Ke;
-}
-/*}}}*/
-/*FUNCTION Seg::CreateKMatrixFreeSurfaceBase {{{*/
-ElementMatrix* Seg::CreateKMatrixFreeSurfaceBase(void){
-
-        /*Intermediaries */
-        int        stabilization;
-        IssmDouble Jdet,D_scalar,dt,h;
-        IssmDouble vx,vel;
-        IssmDouble xyz_list[NUMVERTICES][3];
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = this->NumberofNodes();
-
-        /*Initialize Element matrix and vectors*/
-        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
-        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
-        IssmDouble*    B      = xNew<IssmDouble>(1*numnodes);
-        IssmDouble*    Bprime = xNew<IssmDouble>(1*numnodes);
-
-        /*Retrieve all inputs and parameters*/
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
-        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-        h=this->GetSize();
-
-        /* Start  looping on the number of gaussian points: */
-        GaussSeg *gauss=new GaussSeg(2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-                GetNodalFunctions(basis,gauss);
-
-                vx_input->GetInputValue(&vx,gauss);
-
-                D_scalar=gauss->weight*Jdet;
-
-                TripleMultiply(basis,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        basis,1,numnodes,0,
-                                        &Ke->values[0],1);
-
-                GetNodalFunctions(B,gauss);
-                GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
-
-                D_scalar=dt*gauss->weight*Jdet*vx;
-                TripleMultiply(B,1,numnodes,1,
-                                        &D_scalar,1,1,0,
-                                        Bprime,1,numnodes,0,
-                                        &Ke->values[0],1);
-
-                if(stabilization==2){
-                        /*Streamline upwinding*/
-                        vel=fabs(vx)+1.e-8;
-                        D_scalar=dt*gauss->weight*Jdet*h/(2.*vel)*vx;
-                }
-                else if(stabilization==1){
-                        /*SSA*/
-                        vx_input->GetInputAverage(&vx);
-                        D_scalar=dt*gauss->weight*Jdet*h/2.*fabs(vx);
-                }
-                if(stabilization==1 || stabilization==2){
-                        TripleMultiply(Bprime,1,numnodes,1,
-                                                &D_scalar,1,1,0,
-                                                Bprime,1,numnodes,0,
-                                                &Ke->values[0],1);
-                }
-        }
-
-        /*Clean up and return*/
-        xDelete<IssmDouble>(basis);
-        xDelete<IssmDouble>(B);
-        xDelete<IssmDouble>(Bprime);
-        delete gauss;
-        return Ke;
-}
-/*}}}*/
-/*FUNCTION Seg::CreateMassMatrix {{{*/
-ElementMatrix* Seg::CreateMassMatrix(void){
-
-        /* Intermediaries */
-        IssmDouble  D,Jdet;
-        IssmDouble  xyz_list[NUMVERTICES][3];
-
-        /*Fetch number of nodes and dof for this finite element*/
-        int numnodes = this->NumberofNodes();
-
-        /*Initialize Element matrix and vectors*/
-        ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
-        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-        /*Retrieve all inputs and parameters*/
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-        /* Start looping on the number of gaussian points: */
-        GaussSeg* gauss=new GaussSeg(2);
-        for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                GetNodalFunctions(basis,gauss);
-                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-                D=gauss->weight*Jdet;
-
-                TripleMultiply(basis,1,numnodes,1,
-                                        &D,1,1,0,
-                                        basis,1,numnodes,0,
-                                        &Ke->values[0],1);
-        }
-
-        /*Clean up and return*/
-        delete gauss;
-        xDelete<IssmDouble>(basis);
-        return Ke;
-}
-/*}}}*/
 /*FUNCTION Seg::GetNumberOfNodes;{{{*/
 int Seg::GetNumberOfNodes(void){