Changeset 600

Timestamp:

05/26/09 15:28:42 (16 years ago)

Author:

seroussi

Message:

added PVector in prognostic

Location:

issm/trunk/src/c/objects

Files:

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

issm/trunk/src/c/objects/Tria.cpp

@@ -579 +579 @@
 
         /*recover extra inputs from users, at current convergence iteration: */
-        found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
+        found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
         if(!found)throw ErrorException(__FUNCT__," could not find velocity_average  in inputs!");
 
@@ -728 +728 @@
 }
 
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::CreatePVectorPrognostic"
+void  Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
+
+
+        /* local declarations */
+        int             i,j;
+
+        /* node data: */
+        const int    numgrids=3;
+        const int    NDOF1=1;
+        const int    numdof=NDOF1*numgrids;
+        double       xyz_list[numgrids][3];
+        int          doflist[numdof];
+        int          numberofdofspernode;
+        
+        /* gaussian points: */
+        int     num_gauss,ig;
+        double* first_gauss_area_coord  =  NULL;
+        double* second_gauss_area_coord =  NULL;
+        double* third_gauss_area_coord  =  NULL;
+        double* gauss_weights           =  NULL;
+        double  gauss_weight;
+        double  gauss_l1l2l3[3];
+
+        /* matrix */
+        double pe_g[numgrids]={0.0};
+        double L[numgrids];
+        double Jdettria;
+        
+        /*input parameters for structural analysis (diagnostic): */
+        double  accumulation_list[numgrids]={0.0};
+        double  accumulation_g;
+        double  melting_list[numgrids]={0.0};
+        double  melting_g;
+        double  thickness_list[numgrids]={0.0};
+        double  thickness_g;
+        double  dt;
+        int     dofs[1]={1};
+        int     found=0;
+
+        ParameterInputs* inputs=NULL;
+
+        /*recover pointers: */
+        inputs=(ParameterInputs*)vinputs;
+
+        /*recover extra inputs from users, at current convergence iteration: */
+        found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!");
+        found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!");
+        found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes);
+        if(!found)throw ErrorException(__FUNCT__," could not find thickness in inputs!");
+        found=inputs->Recover("dt",&dt);
+        if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
+
+        /* Get node coordinates and dof list: */
+        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
+        GetDofList(&doflist[0],&numberofdofspernode);
+
+        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+        /* Start  looping on the number of gaussian points: */
+        for (ig=0; ig<num_gauss; ig++){
+                /*Pick up the gaussian point: */
+                gauss_weight=*(gauss_weights+ig);
+                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+                /* Get Jacobian determinant: */
+                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+
+                /*Get L matrix: */
+                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+                /* Get accumulation, melting and thickness at gauss point */
+                GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3);
+                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
+                GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
+                
+                /* Add value into pe_g: */
+                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
+                
+        } // for (ig=0; ig<num_gauss; ig++)
+
+        /*Add pe_g to global matrix Kgg: */
+        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
+
+        cleanup_and_return: 
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&gauss_weights);
+
+}
+
 
 #undef __FUNCT__ 
@@ -1000 +1098 @@
                 
                 CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
+        }
+        else if (analysis_type==PrognosticAnalysisEnum()){
+
+                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
         }
         else{

issm/trunk/src/c/objects/Tria.h

@@ -115 +115 @@
                 void  CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
                 void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
+                void  CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);