Changeset 16936 for issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

Timestamp:

11/25/13 14:38:24 (11 years ago)

Author:

seroussi

Message:

CHG: KMatric coupling HO viscous

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)

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

@@ -828 +828 @@
                 case SSAHOApproximationEnum: 
                         return CreateKMatrixSSAHO(element);
+                case HOFSApproximationEnum: 
+                        return CreateKMatrixHOFS(element);
                 case NoneApproximationEnum:
                         return NULL;
@@ -3065 +3067 @@
 
 /*Coupling (Tiling)*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=CreateKMatrixFS(element);
+        int indices[3]={18,19,20};
+        Ke1->StaticCondensation(3,&indices[0]);
+        int init = element->FiniteElement();
+        element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+        ElementMatrix* Ke2=CreateKMatrixHO(element);
+        element->SetTemporaryElementType(init); // P1 needed for HO
+        ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element);
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
+
+        /*clean-up and return*/
+        delete Ke1;
+        delete Ke2;
+        delete Ke3;
+        return Ke;
+}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
 
@@ -3077 +3098 @@
         delete Ke2;
         delete Ke3;
+        return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingHOFS(Element* element){/*{{{*/
+
+        /*Constants*/
+        int numnodes       = 3*6+1;
+        int numdofp        = 2*6;
+        int numdofs        = 4*6 + 3;
+        int numdoftotal    = (2+4)*6+ 3;
+
+        /*Intermediaries*/
+        int   i,j,init;
+        int*   cs_list   = xNew<int>(6*3+1);
+        int*   cs_list2  = xNew<int>(6*2+1);
+        Node  **node_list = xNew<Node*>(6*3+1);
+
+        /*Some parameters needed*/
+        init = element->FiniteElement();
+
+        /*prepare node list*/
+        for(i=0;i<6+1;i++){
+                node_list[i+6] = element->GetNode(i);
+                cs_list[i+6]   = XYZEnum;
+                cs_list2[i]    = XYZEnum;
+        }
+        for(i=0;i<6;i++){
+                node_list[i]       = element->GetNode(i);
+                node_list[i+2*6+1] = element->GetNode(i+6*1);
+                cs_list[i]         = XYEnum;
+                cs_list[i+2*6+1]   = PressureEnum;
+                cs_list2[i+6+1]    = PressureEnum;
+        }
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=element->NewElementMatrix(FSvelocityEnum);
+        ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
+        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+        delete Ke1; delete Ke2;
+
+        /*Compute HO Matrix with P1 element type\n");*/
+        Ke1=CreateKMatrixFS(element); element->TransformInvStiffnessMatrixCoord(Ke1,node_list,2*6+1,cs_list2);
+        int indices[3]={18,19,20};
+        Ke1->StaticCondensation(3,&indices[0]);
+        element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+        Ke2=CreateKMatrixHO(element); element->TransformInvStiffnessMatrixCoord(Ke2,XYEnum);
+        element->SetTemporaryElementType(init); // P1 needed for HO
+        /*Compute FS Matrix and condense it \n");*/
+
+        for(i=0;i<numdofs;i++) for(j=0;j<6;j++){
+                Ke->values[i*numdoftotal+numdofs+2*j+0]+=Ke1->values[i*numdofs+3*j+0];
+                Ke->values[i*numdoftotal+numdofs+2*j+1]+=Ke1->values[i*numdofs+3*j+1];
+        }
+        for(i=0;i<numdofp;i++) for(j=0;j<6;j++){
+                Ke->values[(i+numdofs)*numdoftotal+3*j+0]+=Ke2->values[i*numdofp+2*j+0];
+                Ke->values[(i+numdofs)*numdoftotal+3*j+1]+=Ke2->values[i*numdofp+2*j+1];
+        }
+
+        /*Transform Coordinate System*/ //Do not transform, already done in the matrices
+        element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+
+        /*clean-up and return*/
+        xDelete<int>(cs_list);
+        xDelete<int>(cs_list2);
+        xDelete<Node*>(node_list);
+        delete Ke1;
+        delete Ke2;
         return Ke;
 }/*}}}*/