Changeset 4942

Timestamp:

08/03/10 09:04:37 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added new functionnalities for ice front

Location:

issm/trunk/src/c/objects

Files:

• Elements/BeamRef.cpp (modified) (1 diff)
• Elements/BeamRef.h (modified) (1 diff)
• Elements/Element.h (modified) (1 diff)
• Elements/Penta.cpp (modified) (3 diffs)
• Elements/Penta.h (modified) (1 diff)
• Elements/Tria.h (modified) (1 diff)
• Loads/Numericalflux.cpp (modified) (2 diffs)

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

@@ -60 +60 @@
 }
 /*}}}*/
+/*FUNCTION BeamRef::GetJacobianDeterminant2d {{{1*/
+void BeamRef::GetJacobianDeterminant2d(double*  Jdet, double* xyz_list,double gauss){
+
+        /*Jdet = 0.5 * length*/
+        *Jdet=1.0/2.0*(sqrt(pow(xyz_list[3*1+0]-xyz_list[3*0+0],2)+pow(xyz_list[3*1+1]-xyz_list[3*0+1],2)));
+        if(*Jdet<0) ISSMERROR(" negative jacobian determinant!");
+
+}
+/*}}}*/
 /*FUNCTION BeamRef::GetNodalFunctions {{{1*/
 void BeamRef::GetNodalFunctions(double* l1l2, double gauss){

issm/trunk/src/c/objects/Elements/BeamRef.h

@@ -24 +24 @@
                 /*Numerics*/
                 void GetJacobianDeterminant(double* Jdet, double* z_list,double gauss);
+                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double gauss);
                 void GetNodalFunctions(double* l1l2, double gauss);
                 void GetParameterValue(double* pvalue, double* value_list,double gauss_coord);

issm/trunk/src/c/objects/Elements/Element.h

@@ -36 +36 @@
                 virtual bool   GetOnBed()=0;
                 virtual void   GetThicknessList(double* thickness_list)=0;
+                virtual void   GetParameterValue(double* pvalue,Node* node,int enumtype)=0;
+                virtual void   GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in)=0;
                 virtual void   GetBedList(double* bed_list)=0;
                 virtual void   Gradj(Vec gradient,int control_type)=0;

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

@@ -4148 +4148 @@
                 delete tria;
         }
-        extern int my_rank;
 
         xfree((void**)&first_gauss_area_coord);
@@ -4236 +4235 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype) {{{1*/
+/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
 void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){
 
@@ -4264 +4263 @@
         return;
 
+}
+/*}}}*/
+/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/
+void Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){
+
+        /*Output*/
+        double value;
+
+        /*Intermediaries*/
+        const int numnodes=6;
+        int       grid1=-1,grid2=-1;
+        int       grid3;
+        int       i;
+        double    gauss_penta[numnodes];
+
+        /*go through 6 nodes (all nodes for penta) and identify 1st and 2nd nodes: */
+        ISSMASSERT(nodes);
+        for(i=0;i<numnodes;i++){
+                if (node1==nodes[i]) grid1=i;
+                if (node2==nodes[i]) grid2=i;
+        }
+
+        /*Reverse grid1 and 2 if necessary*/
+        if (grid1>grid2){
+
+                /*Reverse grid1 and grid2*/
+                grid3=grid1; grid1=grid2; grid2=grid3;
+
+                /*Change segment gauss point (between -1 and +1)*/
+                gauss_seg=-gauss_seg;
+        }
+
+        /*Build Penta Gauss point*/
+        if (grid2<3){
+
+                /*We are on the basal triangle*/
+                gauss_penta[3]=-1.;
+
+                if (grid1==0 && grid2==1){
+                        /*gauss_seg is between 0 and 1*/
+                        gauss_penta[0]=0.5*(1.-gauss_seg);
+                        gauss_penta[1]=1.-0.5*(1.-gauss_seg);
+                        gauss_penta[2]=0.;
+                }
+                else if (grid1==0 && grid2==2){
+                        /*gauss_seg is between 0 and 2*/
+                        gauss_penta[0]=0.5*(1.-gauss_seg);
+                        gauss_penta[1]=0.;
+                        gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+                }
+                else if (grid1==1 && grid2==2){
+                        /*gauss_seg is between 1 and 2*/
+                        gauss_penta[0]=0.;
+                        gauss_penta[1]=0.5*(1.-gauss_seg);
+                        gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+                }
+                else{
+                        ISSMERROR("The 2 nodes provided are either the same or did not match current Penta nodes");
+                }
+        }
+        else if(grid1>2){
+
+                /*We are on the surface triangle*/
+                gauss_penta[3]=+1.;
+
+                if (grid1==3 && grid2==4){
+                        /*gauss_seg is between 0 and 1*/
+                        gauss_penta[0]=0.5*(1.-gauss_seg);
+                        gauss_penta[1]=1.-0.5*(1.-gauss_seg);
+                        gauss_penta[2]=0.;
+                }
+                else if (grid1==3 && grid2==5){
+                        /*gauss_seg is between 0 and 2*/
+                        gauss_penta[0]=0.5*(1.-gauss_seg);
+                        gauss_penta[1]=0.;
+                        gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+                }
+                else if (grid1==4 && grid2==5){
+                        /*gauss_seg is between 1 and 2*/
+                        gauss_penta[0]=0.;
+                        gauss_penta[1]=0.5*(1.-gauss_seg);
+                        gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+                }
+                else{
+                        ISSMERROR("The 2 nodes provided are either the same or did not match current Penta nodes");
+                }
+        }
+        else{
+                ISSMERROR("vertical segments not implemented yet");
+        }
+
+        /*OK, now we can call input method*/
+        input_in->GetParameterValue(&value, &gauss_penta[0]);
+
+        /*Assign output pointers:*/
+        *pvalue=value;
 }
 /*}}}*/

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

@@ -147 +147 @@
                 void      GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
                 void    GetParameterValue(double* pvalue,Node* node,int enumtype);
+                void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
                 void      GetPhi(double* phi, double*  epsilon, double viscosity);
                 void      GetSolutionFromInputsDiagnosticHoriz(Vec solutiong);

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

@@ -142 +142 @@
                 void      GetDofList(int* doflist,int* pnumberofdofs);
                 void      GetDofList1(int* doflist);
-                void      GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
                 void    GetParameterValue(double* pvalue,Node* node,int enumtype);
                 void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);

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

@@ -880 +880 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::GetParameterValue {{{1*/
+/*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype) {{{1*/
 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype){
 
@@ -901 +901 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::GetParameterValue {{{1*/
+/*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/
 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){