Changeset 16936

Timestamp:

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

Author:

seroussi

Message:

CHG: KMatric coupling HO viscous

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (3 diffs)
• analyses/StressbalanceAnalysis.h (modified) (1 diff)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.cpp (modified) (1 diff)

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;
 }/*}}}*/

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

@@ -71 +71 @@
                 void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
                 /*Coupling*/
+                ElementMatrix* CreateKMatrixHOFS(Element* element);
                 ElementMatrix* CreateKMatrixSSAHO(Element* element);
+                ElementMatrix* CreateKMatrixCouplingHOFS(Element* element);
                 ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
                 ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);

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

@@ -296 +296 @@
         for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
 }/*}}}*/
+void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
+        ::TransformInvStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum);
+}/*}}}*/
+void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/
+        ::TransformInvStiffnessMatrixCoord(Ke,nodes_list,numnodes,transformenum_list);
+}/*}}}*/
 void Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/
         ::TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),transformenum);

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -73 +73 @@
                 void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
                 void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
+                void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
+                void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list);
                 void       TransformLoadVectorCoord(ElementVector* pe,int transformenum);
                 void       TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list);

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

@@ -7054 +7054 @@
         /*Compute HO Matrix with P1 element type\n");*/
         this->element_type=P1Enum;
-        Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
+        Ke1=CreateKMatrixStressbalanceHO(); TransformInvStiffnessMatrixCoord(Ke1,XYEnum);
         this->element_type=init;
         /*Compute FS Matrix and condense it \n");*/
-        Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,2*NUMVERTICES+1,cs_list2);
+        Ke2=CreateKMatrixStressbalanceFS(); TransformInvStiffnessMatrixCoord(Ke2,node_list,2*NUMVERTICES+1,cs_list2);
         int indices[3]={18,19,20};
         Ke2->StaticCondensation(3,&indices[0]);