Changeset 5510

Timestamp:

08/23/10 11:01:55 (15 years ago)

Author:

seroussi

Message:

added penalties for Stokes

Location:

issm/trunk/src/c

Files:

• modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (modified) (3 diffs)
• objects/Elements/Penta.cpp (modified) (1 diff)
• objects/Loads/Pengrid.cpp (modified) (2 diffs)
• objects/Loads/Penpair.cpp (modified) (2 diffs)
• objects/Loads/Penpair.h (modified) (1 diff)

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp

@@ -17 +17 @@
         int segment_width;
         int element;
+        int penpair_ids[2];
         
         /*DataSet*/
@@ -81 +82 @@
         for (i=0;i<iomodel->numberofvertices;i++){
                 
-                if((iomodel->my_vertices[i]==1)){
+                if(iomodel->my_vertices[i]==1 && iomodel->singlenodetoelementconnectivity[i]!=0){
 
                         if ((iomodel->gridonbed[i]) && (iomodel->gridonicesheet[i]) && (iomodel->gridonstokes[i])){
@@ -97 +98 @@
         xfree((void**)&iomodel->singlenodetoelementconnectivity);
 
+        /*Create Penpair for penalties: */
+        IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
+        
+        for(i=0;i<iomodel->numpenalties;i++){
+
+                if(iomodel->my_vertices[(int)iomodel->penalties[2*i+0]-1]){
+
+                        /*In debugging mode, check that the second node is in the same cpu*/
+                        ISSMASSERT(iomodel->my_vertices[(int)iomodel->penalties[2*i+1]-1]);
+
+                        /*Get node ids*/
+                        penpair_ids[0]=iomodel->nodecounter+(int)iomodel->penalties[2*i+0];
+                        penpair_ids[1]=iomodel->nodecounter+(int)iomodel->penalties[2*i+1];
+
+                        /*Create Load*/
+                        loads->AddObject(new Penpair(iomodel->loadcounter+count+1,&penpair_ids[0],DiagnosticStokesAnalysisEnum));
+                        count++;
+                }
+        }
+
+        /*free ressources: */
+        xfree((void**)&iomodel->penalties);
+
         cleanup_and_return:
 

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -2158 +2158 @@
                         }
                         if(!iomodel->pressure){
-                                for(i=0;i<6;i++)nodeinputs[i]=i;
+                                for(i=0;i<6;i++)nodeinputs[i]=0;
                                 if(iomodel->qmu_analysis){
                                         this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs));

issm/trunk/src/c/objects/Loads/Pengrid.cpp

@@ -44 +44 @@
         ISSMASSERT(iomodel->singlenodetoelementconnectivity);
         ISSMASSERT(index>=0 && index<iomodel->numberofvertices);
+        ISSMASSERT(id);
 
         /*id: */
@@ -52 +53 @@
         pengrid_node_id=iomodel->nodecounter+index+1;
         pengrid_element_id=iomodel->singlenodetoelementconnectivity[index];
+        ISSMASSERT(pengrid_element_id);
         pengrid_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
 

issm/trunk/src/c/objects/Loads/Penpair.cpp

@@ -209 +209 @@
                 PenaltyCreateKMatrixDiagnosticHoriz(Kgg,kmax);
         }
+        else if (analysis_type==DiagnosticStokesAnalysisEnum){
+                PenaltyCreateKMatrixDiagnosticStokes(Kgg,kmax);
+        }
         else{
                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
@@ -288 +291 @@
 }
 /*}}}1*/
+/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
+void  Penpair::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax){
+        
+        const int numgrids=2;
+        const int NDOF3=4;
+        const int numdof=numgrids*NDOF3;
+        int*      doflist=NULL;
+        int       numberofdofspernode;
+
+        double Ke[8][8]={0.0};
+        double penalty_offset;
+
+        /*pointers: */
+        Node** nodes=NULL;
+
+        /*Get dof list: */
+        GetDofList(&doflist);
+
+        /*recover pointers: */
+        nodes=(Node**)hnodes->deliverp();
+
+        /*recover parameters: */
+        parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
+
+        //Create elementary matrix: add penalty to 
+        Ke[0][0]=kmax*pow((double)10.0,penalty_offset);
+        Ke[0][4]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[4][0]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[4][4]=kmax*pow((double)10.0,penalty_offset);
+
+        Ke[1][1]=kmax*pow((double)10.0,penalty_offset);
+        Ke[1][5]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[5][1]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[5][5]=kmax*pow((double)10.0,penalty_offset);
+        
+        Ke[2][2]=kmax*pow((double)10.0,penalty_offset);
+        Ke[2][6]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[6][2]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[6][6]=kmax*pow((double)10.0,penalty_offset);
+
+        Ke[3][3]=kmax*pow((double)10.0,penalty_offset);
+        Ke[3][7]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[7][3]=-kmax*pow((double)10.0,penalty_offset);
+        Ke[7][7]=kmax*pow((double)10.0,penalty_offset);
+
+        /*Add Ke to global matrix Kgg: */
+        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
+        
+        /*Free ressources:*/
+        xfree((void**)&doflist);
+
+}
+/*}}}1*/
 /*FUNCTION Penpair::GetDofList {{{1*/
 void  Penpair::GetDofList(int** pdoflist){

issm/trunk/src/c/objects/Loads/Penpair.h

@@ -65 +65 @@
                         /*Penpair management: {{{1*/
                 void  PenaltyCreateKMatrixDiagnosticHoriz(Mat Kgg,double kmax);
+                void  PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax);
                 void  GetDofList(int** pdoflist);
                 /*}}}*/