Changeset 5743
- Timestamp:
- 09/10/10 10:20:08 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 25 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Container/Inputs.cpp
r5682 r5743 43 43 44 44 /*Object management*/ 45 /*FUNCTION Inputs::GetParameterValue(double* pvalue, double* gauss,int enum_type){{{1*/46 void Inputs::GetParameterValue(double* pvalue, double* gauss, int enum_type){45 /*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ 46 void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){ 47 47 48 48 vector<Object*>::iterator object; … … 72 72 /*}}}*/ 73 73 /*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ 74 void Inputs::GetParameterValue(double* pvalue,Gauss Tria* gauss, int enum_type){74 void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){ 75 75 76 76 vector<Object*>::iterator object; … … 99 99 } 100 100 /*}}}*/ 101 /*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/102 void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){103 104 vector<Object*>::iterator object;105 Input* input=NULL;106 bool found=false;107 108 /*Go through inputs and check whether any input with the same name is already in: */109 for ( object=objects.begin() ; object < objects.end(); object++ ){110 111 input=(Input*)(*object);112 if (input->EnumType()==enum_type){113 found=true;114 break;115 }116 }117 118 if (!found){119 /*we could not find an input with the correct enum type. No defaults values were provided,120 * error out: */121 ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));122 }123 124 /*Ok, we have an input if we made it here, request the input to return the values: */125 input->GetParameterValue(pvalue,gauss);126 127 }128 /*}}}*/129 /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/130 void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){131 132 vector<Object*>::iterator object;133 Input* input=NULL;134 bool found=false;135 136 /*Go through inputs and check whether any input with the same name is already in: */137 for ( object=objects.begin() ; object < objects.end(); object++ ){138 139 input=(Input*)(*object);140 if (input->EnumType()==enum_type){141 found=true;142 break;143 }144 }145 146 if (!found){147 /*we could not find an input with the correct enum type. Return the default value: */148 *pvalue=defaultvalue;149 }150 else{151 input->GetParameterValue(pvalue,gauss);152 }153 }154 /*}}}*/155 101 /*FUNCTION Inputs::GetParameterValue(bool* pvalue,int enum-type){{{1*/ 156 102 void Inputs::GetParameterValue(bool* pvalue,int enum_type){ … … 263 209 input->GetParameterAverage(pvalue); 264 210 265 }266 /*}}}*/267 /*FUNCTION Inputs::GetParameterDerivativeValue{{{1*/268 void Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){269 270 vector<Object*>::iterator object;271 Input* input=NULL;272 bool found=false;273 274 /*Go through inputs and check whether any input with the same name is already in: */275 for ( object=objects.begin() ; object < objects.end(); object++ ){276 277 input=(Input*)(*object);278 if (input->EnumType()==enum_type){279 found=true;280 break;281 }282 }283 284 if (!found){285 /*we could not find an input with the correct enum type. No defaults values were provided,286 * error out: */287 ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));288 }289 290 /*Ok, we have an input if we made it here, request the input to return the value: */291 input->GetParameterDerivativeValue(derivativevalues,xyz_list,gauss);292 211 } 293 212 /*}}}*/ -
issm/trunk/src/c/Container/Inputs.h
r5682 r5743 48 48 void GetParameterValue(int* pvalue,int enum_type); 49 49 void GetParameterValue(double* pvalue,int enum_type); 50 void GetParameterValue(double* pvalue,double* gauss,int enum_type);51 50 void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type); 52 51 void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type); 53 void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type);55 52 /*}}}*/ 56 53 -
issm/trunk/src/c/objects/Elements/Element.h
r5727 r5743 31 31 virtual void CreatePVector(Vec pg)=0; 32 32 virtual void GetSolutionFromInputs(Vec solution)=0; 33 virtual void GetNodes(void** nodes)=0;34 33 virtual int GetNodeIndex(Node* node)=0; 35 34 virtual bool GetShelf()=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r5736 r5743 884 884 this->results=new Results(); 885 885 886 }887 /*}}}*/888 /*FUNCTION Penta::GetNodes {{{1*/889 void Penta::GetNodes(void** vpnodes){890 891 int i;892 Node** pnodes=NULL;893 894 /*virtual object: */895 pnodes=(Node**)vpnodes;896 897 for(i=0;i<6;i++){898 pnodes[i]=nodes[i];899 }900 886 } 901 887 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Penta.h
r5736 r5743 78 78 void CreatePVector(Vec pg); 79 79 void DeleteResults(void); 80 void GetNodes(void** nodes);81 80 int GetNodeIndex(Node* node); 82 81 bool GetOnBed(); … … 149 148 void GetDofList(int** pdoflist,int approximation_enum=0); 150 149 void GetDofList1(int* doflist); 151 int GetElementType(void); 152 void GetNodalFunctions(double* l1l6, double* gauss_coord); 153 void GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord); 154 void GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord); 155 void GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord); 156 void GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord); 157 void GetNodalFunctionsStokes(double* l1l7, double* gauss_coord); 150 int GetElementType(void); 158 151 void GetParameterListOnVertices(double* pvalue,int enumtype); 159 152 void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r5726 r5743 53 53 /*Reference Element numerics*/ 54 54 /*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/ 55 void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss){55 void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){ 56 56 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 57 57 * For grid i, Bi can be expressed in the actual coordinate system … … 91 91 /*}}}*/ 92 92 /*FUNCTION PentaRef::GetBPattyn {{{1*/ 93 void PentaRef::GetBPattyn(double* B, double* xyz_list, double* gauss){93 void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){ 94 94 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 95 95 * For grid i, Bi can be expressed in the actual coordinate system … … 136 136 /*}}}*/ 137 137 /*FUNCTION PentaRef::GetBprimePattyn {{{1*/ 138 void PentaRef::GetBprimePattyn(double* B, double* xyz_list, double* gauss_coord){138 void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){ 139 139 /*Compute B prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 140 140 * For grid i, Bi can be expressed in the actual coordinate system … … 180 180 /*}}}*/ 181 181 /*FUNCTION PentaRef::GetBStokes {{{1*/ 182 void PentaRef::GetBStokes(double* B, double* xyz_list, double* gauss){182 void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){ 183 183 184 184 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID. … … 251 251 /*}}}*/ 252 252 /*FUNCTION PentaRef::GetBprimeStokes {{{1*/ 253 void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss){253 void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){ 254 254 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 255 255 * For grid i, Bi' can be expressed in the actual coordinate system … … 323 323 /*}}}*/ 324 324 /*FUNCTION PentaRef::GetBArtdiff {{{1*/ 325 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, double* gauss){325 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){ 326 326 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 327 327 * For grid i, Bi' can be expressed in the actual coordinate system … … 353 353 /*}}}*/ 354 354 /*FUNCTION PentaRef::GetBAdvec{{{1*/ 355 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, double* gauss){355 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){ 356 356 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 357 357 * For grid i, Bi' can be expressed in the actual coordinate system … … 385 385 /*}}}*/ 386 386 /*FUNCTION PentaRef::GetBConduct{{{1*/ 387 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, double* gauss){387 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){ 388 388 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 389 389 * For grid i, Bi' can be expressed in the actual coordinate system … … 417 417 /*}}}*/ 418 418 /*FUNCTION PentaRef::GetBVert{{{1*/ 419 void PentaRef::GetBVert(double* B, double* xyz_list, double* gauss){419 void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){ 420 420 /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz]; 421 421 where hi is the interpolation function for grid i.*/ … … 437 437 /*}}}*/ 438 438 /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/ 439 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, double* gauss){439 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){ 440 440 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. 441 441 * For grid i, Bi' can be expressed in the actual coordinate system … … 469 469 /*}}}*/ 470 470 /*FUNCTION PentaRef::GetBprimeVert{{{1*/ 471 void PentaRef::GetBprimeVert(double* B, double* xyz_list, double* gauss){471 void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){ 472 472 /* Compute Bprime matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/ 473 473 … … 477 477 /*}}}*/ 478 478 /*FUNCTION PentaRef::GetLStokes {{{1*/ 479 void PentaRef::GetLStokes(double* LStokes, double* gauss){479 void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){ 480 480 /* 481 481 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. … … 507 507 508 508 /*Get l1l2l3 in actual coordinate system: */ 509 l1l2l3[0]=gauss [0];510 l1l2l3[1]=gauss [1];511 l1l2l3[2]=gauss [2];509 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; 510 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; 511 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; 512 512 513 513 /*Build LStokes: */ … … 574 574 /*}}}*/ 575 575 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/ 576 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss){576 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){ 577 577 578 578 /* … … 604 604 605 605 /*Get l1l2l3 in actual coordinate system: */ 606 l1l2l3[0]=gauss [0];607 l1l2l3[1]=gauss [1];608 l1l2l3[2]=gauss [2];606 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0; 607 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0; 608 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0; 609 609 610 610 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss); … … 673 673 /*}}}*/ 674 674 /*FUNCTION PentaRef::GetJacobian {{{1*/ 675 void PentaRef::GetJacobian(double* J, double* xyz_list,double* gauss){676 677 const int NDOF3=3;678 int i,j;679 680 /*The Jacobian is constant over the element, discard the gaussian points.681 * J is assumed to have been allocated of size NDOF2xNDOF2.*/682 683 double A1,A2,A3; //area coordinates684 double xi,eta,zi; //parametric coordinates685 686 double x1,x2,x3,x4,x5,x6;687 double y1,y2,y3,y4,y5,y6;688 double z1,z2,z3,z4,z5,z6;689 690 /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */691 A1=gauss[0];692 A2=gauss[1];693 A3=gauss[2];694 695 xi=A2-A1;696 eta=SQRT3*A3;697 zi=gauss[3];698 699 x1=*(xyz_list+3*0+0);700 x2=*(xyz_list+3*1+0);701 x3=*(xyz_list+3*2+0);702 x4=*(xyz_list+3*3+0);703 x5=*(xyz_list+3*4+0);704 x6=*(xyz_list+3*5+0);705 706 y1=*(xyz_list+3*0+1);707 y2=*(xyz_list+3*1+1);708 y3=*(xyz_list+3*2+1);709 y4=*(xyz_list+3*3+1);710 y5=*(xyz_list+3*4+1);711 y6=*(xyz_list+3*5+1);712 713 z1=*(xyz_list+3*0+2);714 z2=*(xyz_list+3*1+2);715 z3=*(xyz_list+3*2+2);716 z4=*(xyz_list+3*3+2);717 z5=*(xyz_list+3*4+2);718 z6=*(xyz_list+3*5+2);719 720 *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);721 *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);722 *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);723 724 *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);725 *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);726 *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);727 728 *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);729 *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);730 *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);731 732 }733 /*}}}*/734 /*FUNCTION PentaRef::GetJacobianDeterminant {{{1*/735 void PentaRef::GetJacobianDeterminant(double* Jdet, double* xyz_list,double* gauss){736 /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take737 * the determinant of it: */738 double J[3][3];739 740 /*Get Jacobian*/741 GetJacobian(&J[0][0],xyz_list,gauss);742 743 /*Get Determinant*/744 Matrix3x3Determinant(Jdet,&J[0][0]);745 if(*Jdet<0) ISSMERROR("negative jacobian determinant!");746 747 }748 /*}}}*/749 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/750 void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){751 752 /*Jacobian*/753 double J[3][3];754 755 /*Call Jacobian routine to get the jacobian:*/756 GetJacobian(&J[0][0], xyz_list, gauss);757 758 /*Invert Jacobian matrix: */759 Matrix3x3Invert(Jinv,&J[0][0]);760 }761 /*}}}*/762 /*FUNCTION PentaRef::GetNodalFunctionsMINI{{{1*/763 void PentaRef::GetNodalFunctionsMINI(double* l1l7, double* gauss){764 /*This routine returns the values of the nodal functions at the gaussian point.*/765 766 l1l7[0]=gauss[0]*(1.0-gauss[3])/2.0;767 l1l7[1]=gauss[1]*(1.0-gauss[3])/2.0;768 l1l7[2]=gauss[2]*(1.0-gauss[3])/2.0;769 l1l7[3]=gauss[0]*(1.0+gauss[3])/2.0;770 l1l7[4]=gauss[1]*(1.0+gauss[3])/2.0;771 l1l7[5]=gauss[2]*(1.0+gauss[3])/2.0;772 l1l7[6]=27*gauss[0]*gauss[1]*gauss[2]*(1.0+gauss[3])*(1.0-gauss[3]);773 774 }775 /*}}}*/776 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{1*/777 void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss){778 779 /*This routine returns the values of the nodal functions derivatives (with respect to the780 * actual coordinate system): */781 782 int i;783 const int numgrids = 7;784 double dh1dh7_ref[3][numgrids];785 double Jinv[3][3];786 787 /*Get derivative values with respect to parametric coordinate system: */788 GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss);789 790 /*Get Jacobian invert: */791 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);792 793 /*Build dh1dh6:794 *795 * [dhi/dx]= Jinv'*[dhi/dr]796 * [dhi/dy] [dhi/ds]797 * [dhi/dz] [dhi/dzeta]798 */799 800 for (i=0;i<numgrids;i++){801 *(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];802 *(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];803 *(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];804 }805 806 }807 /*}}}*/808 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{1*/809 void PentaRef::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss){810 811 /*This routine returns the values of the nodal functions derivatives (with respect to the812 * natural coordinate system) at the gaussian point. */813 814 int numgrids=7; //six plus bubble grids815 816 double r=gauss[1]-gauss[0];817 double s=-3.0/SQRT3*(gauss[0]+gauss[1]-2.0/3.0);818 double zeta=gauss[3];819 820 /*First nodal function: */821 *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;822 *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;823 *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);824 825 /*Second nodal function: */826 *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;827 *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;828 *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);829 830 /*Third nodal function: */831 *(dl1dl7+numgrids*0+2)=0;832 *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;833 *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);834 835 /*Fourth nodal function: */836 *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;837 *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;838 *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);839 840 /*Fith nodal function: */841 *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;842 *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;843 *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);844 845 /*Sixth nodal function: */846 *(dl1dl7+numgrids*0+5)=0;847 *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;848 *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);849 850 /*Seventh nodal function: */851 *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);852 *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));853 *(dl1dl7+numgrids*2+6)=27*gauss[0]*gauss[1]*gauss[2]*(-2.0*zeta);854 855 }856 /*}}}*/857 /*FUNCTION PentaRef::GetNodalFunctionsP1 {{{1*/858 void PentaRef::GetNodalFunctionsP1(double* l1l6, double* gauss){859 /*This routine returns the values of the nodal functions at the gaussian point.*/860 861 l1l6[0]=gauss[0]*(1-gauss[3])/2.0;862 l1l6[1]=gauss[1]*(1-gauss[3])/2.0;863 l1l6[2]=gauss[2]*(1-gauss[3])/2.0;864 l1l6[3]=gauss[0]*(1+gauss[3])/2.0;865 l1l6[4]=gauss[1]*(1+gauss[3])/2.0;866 l1l6[5]=gauss[2]*(1+gauss[3])/2.0;867 868 }869 /*}}}*/870 /*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{1*/871 void PentaRef::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss){872 873 /*This routine returns the values of the nodal functions derivatives (with respect to the874 * actual coordinate system): */875 int i;876 const int NDOF3 = 3;877 const int numgrids = 6;878 double dh1dh6_ref[NDOF3][numgrids];879 double Jinv[NDOF3][NDOF3];880 881 /*Get derivative values with respect to parametric coordinate system: */882 GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss);883 884 /*Get Jacobian invert: */885 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);886 887 /*Build dh1dh3:888 *889 * [dhi/dx]= Jinv*[dhi/dr]890 * [dhi/dy] [dhi/ds]891 * [dhi/dz] [dhi/dn]892 */893 894 for (i=0;i<numgrids;i++){895 *(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];896 *(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];897 *(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];898 }899 900 }901 /*}}}*/902 /*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{1*/903 void PentaRef::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss){904 905 /*This routine returns the values of the nodal functions derivatives (with respect to the906 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */907 908 const int numgrids=6;909 double A1,A2,A3,z;910 911 A1=gauss[0]; ISSMASSERT(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);912 A2=gauss[1]; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);913 A3=gauss[2]; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value In term of xi and eta: A3=y/SQRT3;914 z=gauss[3]; ISSMASSERT(z>=-1 && z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2915 916 /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/917 *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;918 *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;919 *(dl1dl6+numgrids*2+0)=-0.5*A1;920 921 /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/922 *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;923 *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;924 *(dl1dl6+numgrids*2+1)=-0.5*A2;925 926 /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/927 *(dl1dl6+numgrids*0+2)=0.0;928 *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;929 *(dl1dl6+numgrids*2+2)=-0.5*A3;930 931 /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/932 *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;933 *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;934 *(dl1dl6+numgrids*2+3)=0.5*A1;935 936 /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/937 *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;938 *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;939 *(dl1dl6+numgrids*2+4)=0.5*A2;940 941 /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/942 *(dl1dl6+numgrids*0+5)=0.0;943 *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;944 *(dl1dl6+numgrids*2+5)=0.5*A3;945 }946 /*}}}*/947 /*FUNCTION PentaRef::GetParameterValue{{{1*/948 void PentaRef::GetParameterValue(double* pvalue,double* plist,double* gauss){949 /*P1 interpolation on Gauss point*/950 951 /*intermediary*/952 double l1l6[6];953 954 /*nodal functions: */955 GetNodalFunctionsP1(&l1l6[0],gauss);956 957 /*Assign output pointers:*/958 *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5];959 960 }961 /*}}}*/962 /*FUNCTION PentaRef::GetParameterDerivativeValue{{{1*/963 void PentaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss){964 /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:965 * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;966 * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;967 * dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;968 *969 * p is a vector of size 3x1 already allocated.970 */971 double dh1dh6[3][6];972 973 /*Get nodal funnctions derivatives in actual coordinate system: */974 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);975 976 /*Assign output*/977 p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5];978 p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5];979 p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5];980 981 }982 /*}}}*/983 984 /*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/985 void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){986 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.987 * For grid i, Bi can be expressed in the actual coordinate system988 * by:989 * Bi=[ dh/dx 0 ]990 * [ 0 dh/dy ]991 * [ 1/2*dh/dy 1/2*dh/dx ]992 * where h is the interpolation function for grid i.993 *994 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)995 */996 997 int i;998 const int numgrids=6;999 const int NDOF3=3;1000 const int NDOF2=2;1001 1002 double dh1dh6[NDOF3][numgrids];1003 1004 /*Get dh1dh6 in actual coordinate system: */1005 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);1006 1007 /*Build B: */1008 for (i=0;i<numgrids;i++){1009 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];1010 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;1011 1012 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;1013 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];1014 1015 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];1016 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];1017 1018 }1019 1020 }1021 /*}}}*/1022 /*FUNCTION PentaRef::GetBPattyn {{{1*/1023 void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){1024 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.1025 * For grid i, Bi can be expressed in the actual coordinate system1026 * by:1027 * Bi=[ dh/dx 0 ]1028 * [ 0 dh/dy ]1029 * [ 1/2*dh/dy 1/2*dh/dx ]1030 * [ 1/2*dh/dz 0 ]1031 * [ 0 1/2*dh/dz ]1032 * where h is the interpolation function for grid i.1033 *1034 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)1035 */1036 1037 int i;1038 const int numgrids=6;1039 const int NDOF3=3;1040 const int NDOF2=2;1041 1042 double dh1dh6[NDOF3][numgrids];1043 1044 /*Get dh1dh6 in actual coordinate system: */1045 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);1046 1047 /*Build B: */1048 for (i=0;i<numgrids;i++){1049 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];1050 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;1051 1052 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;1053 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];1054 1055 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];1056 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];1057 1058 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];1059 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;1060 1061 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;1062 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];1063 }1064 1065 }1066 /*}}}*/1067 /*FUNCTION PentaRef::GetBprimePattyn {{{1*/1068 void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){1069 /*Compute B prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.1070 * For grid i, Bi can be expressed in the actual coordinate system1071 * by:1072 * Bi=[ 2*dh/dx dh/dy ]1073 * [ dh/dx 2*dh/dy ]1074 * [ dh/dy dh/dx ]1075 * [ dh/dz 0 ]1076 * [ 0 dh/dz ]1077 * where h is the interpolation function for grid i.1078 *1079 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)1080 */1081 1082 int i;1083 const int NDOF3=3;1084 const int NDOF2=2;1085 const int numgrids=6;1086 1087 double dh1dh6[NDOF3][numgrids];1088 1089 /*Get dh1dh6 in actual coordinate system: */1090 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);1091 1092 /*Build BPrime: */1093 for (i=0;i<numgrids;i++){1094 *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i];1095 *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];1096 1097 *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];1098 *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];1099 1100 *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i];1101 *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i];1102 1103 *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i];1104 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;1105 1106 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;1107 *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i];1108 }1109 }1110 /*}}}*/1111 /*FUNCTION PentaRef::GetBStokes {{{1*/1112 void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){1113 1114 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.1115 * For grid i, Bi can be expressed in the actual coordinate system1116 * by: Bi=[ dh/dx 0 0 0 ]1117 * [ 0 dh/dy 0 0 ]1118 * [ 0 0 dh/dy 0 ]1119 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ]1120 * [ 1/2*dh/dz 0 1/2*dh/dx 0 ]1121 * [ 0 1/2*dh/dz 1/2*dh/dy 0 ]1122 * [ 0 0 0 h ]1123 * [ dh/dx dh/dy dh/dz 0 ]1124 * where h is the interpolation function for grid i.1125 * Same thing for Bb except the last column that does not exist.1126 */1127 1128 int i;1129 const int calculationdof=3;1130 const int numgrids=6;1131 int DOFPERGRID=4;1132 1133 double dh1dh7[calculationdof][numgrids+1];1134 double l1l6[numgrids];1135 1136 1137 /*Get dh1dh7 in actual coordinate system: */1138 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);1139 GetNodalFunctionsP1(l1l6, gauss);1140 1141 /*Build B: */1142 for (i=0;i<numgrids+1;i++){1143 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];1144 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;1145 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;1146 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;1147 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];1148 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;1149 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;1150 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;1151 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];1152 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];1153 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];1154 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;1155 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];1156 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;1157 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];1158 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;1159 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];1160 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];1161 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;1162 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;1163 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;1164 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];1165 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];1166 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];1167 }1168 1169 for (i=0;i<numgrids;i++){ //last column not for the bubble function1170 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;1171 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;1172 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;1173 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;1174 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;1175 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;1176 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];1177 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;1178 }1179 1180 }1181 /*}}}*/1182 /*FUNCTION PentaRef::GetBprimeStokes {{{1*/1183 void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){1184 /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.1185 * For grid i, Bi' can be expressed in the actual coordinate system1186 * by:1187 * Bi'=[ dh/dx 0 0 0]1188 * [ 0 dh/dy 0 0]1189 * [ 0 0 dh/dz 0]1190 * [ dh/dy dh/dx 0 0]1191 * [ dh/dz 0 dh/dx 0]1192 * [ 0 dh/dz dh/dy 0]1193 * [ dh/dx dh/dy dh/dz 0]1194 * [ 0 0 0 h]1195 * where h is the interpolation function for grid i.1196 *1197 * Same thing for the bubble fonction except that there is no fourth column1198 */1199 1200 int i;1201 const int calculationdof=3;1202 const int numgrids=6;1203 int DOFPERGRID=4;1204 1205 double dh1dh7[calculationdof][numgrids+1];1206 double l1l6[numgrids];1207 1208 /*Get dh1dh7 in actual coordinate system: */1209 GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);1210 1211 GetNodalFunctionsP1(l1l6, gauss);1212 1213 /*B_primeuild B_prime: */1214 for (i=0;i<numgrids+1;i++){1215 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];1216 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;1217 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;1218 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;1219 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];1220 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;1221 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;1222 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;1223 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];1224 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i];1225 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i];1226 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;1227 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];1228 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;1229 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];1230 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;1231 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];1232 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];1233 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];1234 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];1235 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];1236 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;1237 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;1238 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;1239 }1240 1241 for (i=0;i<numgrids;i++){ //last column not for the bubble function1242 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;1243 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;1244 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;1245 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;1246 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;1247 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;1248 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;1249 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];1250 }1251 1252 }1253 /*}}}*/1254 /*FUNCTION PentaRef::GetBArtdiff {{{1*/1255 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){1256 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.1257 * For grid i, Bi' can be expressed in the actual coordinate system1258 * by:1259 * Bi_artdiff=[ dh/dx ]1260 * [ dh/dy ]1261 * where h is the interpolation function for grid i.1262 *1263 * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)1264 */1265 1266 int i;1267 const int calculationdof=3;1268 const int numgrids=6;1269 int DOFPERGRID=1;1270 1271 /*Same thing in the actual coordinate system: */1272 double dh1dh6[calculationdof][numgrids];1273 1274 /*Get dh1dh6 in actual coordinates system : */1275 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);1276 1277 /*Build B': */1278 for (i=0;i<numgrids;i++){1279 *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];1280 *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];1281 }1282 }1283 /*}}}*/1284 /*FUNCTION PentaRef::GetBAdvec{{{1*/1285 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){1286 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.1287 * For grid i, Bi' can be expressed in the actual coordinate system1288 * by:1289 * Bi_advec =[ h ]1290 * [ h ]1291 * [ h ]1292 * where h is the interpolation function for grid i.1293 *1294 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)1295 */1296 1297 int i;1298 int calculationdof=3;1299 int numgrids=6;1300 int DOFPERGRID=1;1301 1302 /*Same thing in the actual coordinate system: */1303 double l1l6[6];1304 1305 /*Get dh1dh2dh3 in actual coordinates system : */1306 GetNodalFunctionsP1(l1l6, gauss);1307 1308 /*Build B': */1309 for (i=0;i<numgrids;i++){1310 *(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i];1311 *(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i];1312 *(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i];1313 }1314 }1315 /*}}}*/1316 /*FUNCTION PentaRef::GetBConduct{{{1*/1317 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){1318 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.1319 * For grid i, Bi' can be expressed in the actual coordinate system1320 * by:1321 * Bi_conduct=[ dh/dx ]1322 * [ dh/dy ]1323 * [ dh/dz ]1324 * where h is the interpolation function for grid i.1325 *1326 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)1327 */1328 1329 int i;1330 const int calculationdof=3;1331 const int numgrids=6;1332 int DOFPERGRID=1;1333 1334 /*Same thing in the actual coordinate system: */1335 double dh1dh6[calculationdof][numgrids];1336 1337 /*Get dh1dh2dh3 in actual coordinates system : */1338 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);1339 1340 /*Build B': */1341 for (i=0;i<numgrids;i++){1342 *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];1343 *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];1344 *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];1345 }1346 }1347 /*}}}*/1348 /*FUNCTION PentaRef::GetBVert{{{1*/1349 void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){1350 /* Compute B matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];1351 where hi is the interpolation function for grid i.*/1352 1353 int i;1354 const int NDOF3=3;1355 const int numgrids=6;1356 double dh1dh6[NDOF3][numgrids];1357 1358 /*Get dh1dh6 in actual coordinate system: */1359 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);1360 1361 /*Build B: */1362 for (i=0;i<numgrids;i++){1363 B[i]=dh1dh6[2][i];1364 }1365 1366 }1367 /*}}}*/1368 /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/1369 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){1370 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.1371 * For grid i, Bi' can be expressed in the actual coordinate system1372 * by:1373 * Biprime_advec=[ dh/dx ]1374 * [ dh/dy ]1375 * [ dh/dz ]1376 * where h is the interpolation function for grid i.1377 *1378 * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)1379 */1380 1381 int i;1382 const int calculationdof=3;1383 const int numgrids=6;1384 int DOFPERGRID=1;1385 1386 /*Same thing in the actual coordinate system: */1387 double dh1dh6[calculationdof][numgrids];1388 1389 /*Get dh1dh2dh3 in actual coordinates system : */1390 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);1391 1392 /*Build B': */1393 for (i=0;i<numgrids;i++){1394 *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];1395 *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];1396 *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];1397 }1398 }1399 /*}}}*/1400 /*FUNCTION PentaRef::GetBprimeVert{{{1*/1401 void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){1402 /* Compute Bprime matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/1403 1404 GetNodalFunctionsP1(B, gauss);1405 1406 }1407 /*}}}*/1408 /*FUNCTION PentaRef::GetLStokes {{{1*/1409 void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){1410 /*1411 * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.1412 * For grid i, Li can be expressed in the actual coordinate system1413 * by:1414 * Li=[ h 0 0 0]1415 * [ 0 h 0 0]1416 * [ 0 0 h 0]1417 * [ 0 0 h 0]1418 * [ h 0 0 0]1419 * [ 0 h 0 0]1420 * [ h 0 0 0]1421 * [ 0 h 0 0]1422 * [ 0 0 h 0]1423 * [ 0 0 h 0]1424 * [ 0 0 h 0]1425 * [ h 0 0 0]1426 * [ 0 h 0 0]1427 * [ 0 0 h 0]1428 * where h is the interpolation function for grid i.1429 */1430 1431 int i;1432 const int numgrids2d=3;1433 int num_dof=4;1434 1435 double l1l2l3[numgrids2d];1436 1437 1438 /*Get l1l2l3 in actual coordinate system: */1439 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;1440 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;1441 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;1442 1443 /*Build LStokes: */1444 for (i=0;i<3;i++){1445 *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];1446 *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;1447 *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;1448 *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;1449 *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;1450 *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];1451 *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;1452 *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;1453 *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;1454 *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;1455 *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];1456 *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;1457 *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;1458 *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;1459 *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];1460 *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;1461 *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];1462 *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;1463 *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;1464 *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;1465 *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;1466 *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];1467 *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;1468 *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;1469 *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];1470 *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;1471 *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;1472 *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;1473 *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;1474 *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];1475 *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;1476 *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;1477 *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;1478 *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;1479 *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];1480 *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;1481 *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;1482 *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;1483 *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];1484 *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;1485 *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;1486 *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;1487 *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];1488 *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;1489 *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];1490 *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;1491 *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;1492 *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;1493 *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;1494 *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];1495 *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;1496 *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;1497 *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;1498 *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;1499 *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];1500 *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;1501 1502 }1503 }1504 /*}}}*/1505 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/1506 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){1507 1508 /*1509 * Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.1510 * For grid i, Lpi can be expressed in the actual coordinate system1511 * by:1512 * Lpi=[ h 0 0 0]1513 * [ 0 h 0 0]1514 * [ h 0 0 0]1515 * [ 0 h 0 0]1516 * [ 0 0 h 0]1517 * [ 0 0 h 0]1518 * [ 0 0 dh/dz 0]1519 * [ 0 0 dh/dz 0]1520 * [ 0 0 dh/dz 0]1521 * [dh/dz 0 dh/dx 0]1522 * [ 0 dh/dz dh/dy 0]1523 * [ 0 0 0 h]1524 * [ 0 0 0 h]1525 * [ 0 0 0 h]1526 * where h is the interpolation function for grid i.1527 */1528 int i;1529 const int numgrids2d=3;1530 int num_dof=4;1531 1532 double l1l2l3[numgrids2d];1533 double dh1dh6[3][6];1534 1535 /*Get l1l2l3 in actual coordinate system: */1536 l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;1537 l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;1538 l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;1539 1540 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);1541 1542 /*Build LprimeStokes: */1543 for (i=0;i<3;i++){1544 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];1545 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;1546 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;1547 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;1548 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0;1549 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];1550 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;1551 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;1552 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i];1553 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;1554 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0;1555 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;1556 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0;1557 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i];1558 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0;1559 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;1560 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0;1561 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;1562 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i];1563 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;1564 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0;1565 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0;1566 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i];1567 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;1568 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;1569 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;1570 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];1571 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;1572 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;1573 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;1574 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];1575 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;1576 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;1577 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;1578 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];1579 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;1580 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];1581 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;1582 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];1583 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;1584 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;1585 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];1586 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];1587 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;1588 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;1589 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;1590 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;1591 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i];1592 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0;1593 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0;1594 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;1595 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i];1596 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0;1597 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;1598 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;1599 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];1600 1601 }1602 }1603 /*}}}*/1604 /*FUNCTION PentaRef::GetJacobian {{{1*/1605 675 void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){ 1606 676 -
issm/trunk/src/c/objects/Elements/PentaRef.h
r5726 r5743 23 23 24 24 /*Numerics*/ 25 void GetNodalFunctionsP1(double* l1l6, double* gauss);26 void GetNodalFunctionsMINI(double* l1l7, double* gauss);27 void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss);28 void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss);29 void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss);30 void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss);31 void GetJacobian(double* J, double* xyz_list,double* gauss);32 void GetJacobianDeterminant(double* Jdet, double* xyz_list,double* gauss);33 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss);34 void GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss);35 void GetBPattyn(double* B, double* xyz_list, double* gauss);36 void GetBprimePattyn(double* B, double* xyz_list, double* gauss);37 void GetBStokes(double* B, double* xyz_list, double* gauss);38 void GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss);39 void GetBprimeVert(double* B, double* xyz_list, double* gauss);40 void GetBAdvec(double* B_advec, double* xyz_list, double* gauss);41 void GetBArtdiff(double* B_artdiff, double* xyz_list, double* gauss);42 void GetBConduct(double* B_conduct, double* xyz_list, double* gauss);43 void GetBVert(double* B, double* xyz_list, double* gauss);44 void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, double* gauss);45 void GetLStokes(double* LStokes, double* gauss);46 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss);47 void GetParameterValue(double* pvalue,double* plist,double* gauss);48 void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};49 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, double* gauss);50 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};51 52 25 void GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss); 53 26 void GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss); … … 75 48 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss); 76 49 void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss); 77 //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};50 void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; 78 51 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss); 79 //void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};52 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; 80 53 81 54 }; -
issm/trunk/src/c/objects/Elements/Tria.cpp
r5741 r5743 807 807 } 808 808 /*}}}*/ 809 /*FUNCTION Tria::GetNodes {{{1*/810 void Tria::GetNodes(void** vpnodes){811 int i;812 Node** pnodes=NULL;813 814 /*recover nodes: */815 pnodes=(Node**)vpnodes;816 817 for(i=0;i<3;i++){818 pnodes[i]=nodes[i];819 }820 }821 /*}}}*/822 809 /*FUNCTION Tria::GetNodeIndex {{{1*/ 823 810 int Tria::GetNodeIndex(Node* node){ … … 5502 5489 } 5503 5490 /*}}}*/ 5504 /*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype) {{{1*/5505 void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype){5506 5507 /*Output*/5508 double value;5509 5510 /*Intermediaries*/5511 const int numnodes=3;5512 int grid1=-1,grid2=-1;5513 int grid3;5514 int i;5515 double gauss_tria[numnodes];5516 5517 /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */5518 ISSMASSERT(nodes);5519 for(i=0;i<numnodes;i++){5520 if (node1==nodes[i]) grid1=i;5521 if (node2==nodes[i]) grid2=i;5522 }5523 5524 /*Reverse grid1 and 2 if necessary*/5525 if (grid1>grid2){5526 5527 /*Reverse grid1 and grid2*/5528 grid3=grid1; grid1=grid2; grid2=grid3;5529 5530 /*Change segment gauss point (between -1 and +1)*/5531 gauss_seg=-gauss_seg;5532 }5533 5534 /*Build Triangle Gauss point*/5535 if (grid1==0 && grid2==1){5536 /*gauss_seg is between 0 and 1*/5537 gauss_tria[0]=0.5*(1.-gauss_seg);5538 gauss_tria[1]=1.-0.5*(1.-gauss_seg);5539 gauss_tria[2]=0.;5540 }5541 else if (grid1==0 && grid2==2){5542 /*gauss_seg is between 0 and 2*/5543 gauss_tria[0]=0.5*(1.-gauss_seg);5544 gauss_tria[1]=0.;5545 gauss_tria[2]=1.-0.5*(1.-gauss_seg);5546 }5547 else if (grid1==1 && grid2==2){5548 /*gauss_seg is between 1 and 2*/5549 gauss_tria[0]=0.;5550 gauss_tria[1]=0.5*(1.-gauss_seg);5551 gauss_tria[2]=1.-0.5*(1.-gauss_seg);5552 }5553 else{5554 ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes");5555 }5556 5557 /*OK, now we can call input method*/5558 this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype);5559 5560 /*Assign output pointers:*/5561 *pvalue=value;5562 }5563 /*}}}*/5564 /*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/5565 void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){5566 5567 /*Output*/5568 double value;5569 5570 /*Intermediaries*/5571 const int numnodes=3;5572 int grid1=-1,grid2=-1;5573 int grid3;5574 int i;5575 double gauss_tria[numnodes];5576 5577 /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */5578 ISSMASSERT(nodes);5579 for(i=0;i<numnodes;i++){5580 if (node1==nodes[i]) grid1=i;5581 if (node2==nodes[i]) grid2=i;5582 }5583 5584 /*Reverse grid1 and 2 if necessary*/5585 if (grid1>grid2){5586 5587 /*Reverse grid1 and grid2*/5588 grid3=grid1; grid1=grid2; grid2=grid3;5589 5590 /*Change segment gauss point (between -1 and +1)*/5591 gauss_seg=-gauss_seg;5592 }5593 5594 /*Build Triangle Gauss point*/5595 if (grid1==0 && grid2==1){5596 /*gauss_seg is between 0 and 1*/5597 gauss_tria[0]=0.5*(1.-gauss_seg);5598 gauss_tria[1]=1.-0.5*(1.-gauss_seg);5599 gauss_tria[2]=0.;5600 }5601 else if (grid1==0 && grid2==2){5602 /*gauss_seg is between 0 and 2*/5603 gauss_tria[0]=0.5*(1.-gauss_seg);5604 gauss_tria[1]=0.;5605 gauss_tria[2]=1.-0.5*(1.-gauss_seg);5606 }5607 else if (grid1==1 && grid2==2){5608 /*gauss_seg is between 1 and 2*/5609 gauss_tria[0]=0.;5610 gauss_tria[1]=0.5*(1.-gauss_seg);5611 gauss_tria[2]=1.-0.5*(1.-gauss_seg);5612 }5613 else{5614 ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes");5615 }5616 5617 /*OK, now we can call input method*/5618 input_in->GetParameterValue(&value, &gauss_tria[0]);5619 5620 /*Assign output pointers:*/5621 *pvalue=value;5622 }5623 /*}}}*/5624 5491 /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/ 5625 5492 void Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ -
issm/trunk/src/c/objects/Elements/Tria.h
r5734 r5743 74 74 void CreateKMatrix(Mat Kgg); 75 75 void CreatePVector(Vec pg); 76 void GetNodes(void** nodes);77 76 int GetNodeIndex(Node* node); 78 77 bool GetOnBed(); … … 153 152 void GetParameterListOnVertices(double* pvalue,int enumtype); 154 153 void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); 155 void GetParameterValue(double* pvalue,Node* node,int enumtype); 156 void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype); 157 void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in); 154 void GetParameterValue(double* pvalue,Node* node,int enumtype); 158 155 void GetSolutionFromInputsDiagnosticHoriz(Vec solution); 159 156 void GetSolutionFromInputsDiagnosticHutter(Vec solution); -
issm/trunk/src/c/objects/Elements/TriaRef.cpp
r5738 r5743 52 52 53 53 /*Reference Element numerics*/ 54 /*FUNCTION TriaRef::GetBMacAyeal {{{1*/55 void TriaRef::GetBMacAyeal(double* B, double* xyz_list, double* gauss){56 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.57 * For grid i, Bi can be expressed in the actual coordinate system58 * by:59 * Bi=[ dh/dx 0 ]60 * [ 0 dh/dy ]61 * [ 1/2*dh/dy 1/2*dh/dx ]62 * where h is the interpolation function for grid i.63 *64 * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)65 */66 67 int i;68 const int NDOF2=2;69 const int numgrids=3;70 71 double dh1dh3[NDOF2][numgrids];72 73 74 /*Get dh1dh2dh3 in actual coordinate system: */75 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);76 77 /*Build B: */78 for (i=0;i<numgrids;i++){79 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];80 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;81 *(B+NDOF2*numgrids*1+NDOF2*i)=0;82 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];83 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];84 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];85 }86 }87 /*}}}*/88 54 /*FUNCTION TriaRef::GetBMacAyeal {{{1*/ 89 55 void TriaRef::GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss){ … … 161 127 /*}}}*/ 162 128 /*FUNCTION TriaRef::GetBPrognostic{{{1*/ 163 void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, double* gauss){164 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.165 * For grid i, Bi can be expressed in the actual coordinate system166 * by:167 * Bi=[ h ]168 * [ h ]169 * where h is the interpolation function for grid i.170 *171 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)172 */173 174 int i;175 const int NDOF1=1;176 const int numgrids=3;177 178 double l1l2l3[numgrids];179 180 181 /*Get dh1dh2dh3 in actual coordinate system: */182 GetNodalFunctions(&l1l2l3[0],gauss);183 184 /*Build B_prog: */185 for (i=0;i<numgrids;i++){186 *(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];187 *(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];188 }189 }190 /*}}}*/191 /*FUNCTION TriaRef::GetBPrognostic{{{1*/192 129 void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss){ 193 130 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. … … 219 156 /*}}}*/ 220 157 /*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/ 221 void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss){158 void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){ 222 159 223 160 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. … … 253 190 } 254 191 /*}}}*/ 255 /*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/256 void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){257 258 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.259 * For grid i, Bi' can be expressed in the actual coordinate system260 * by:261 * Bi_prime=[ 2*dh/dx dh/dy ]262 * [ dh/dx 2*dh/dy ]263 * [ dh/dy dh/dx ]264 * where h is the interpolation function for grid i.265 *266 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)267 */268 269 int i;270 const int NDOF2=2;271 const int numgrids=3;272 273 /*Same thing in the actual coordinate system: */274 double dh1dh3[NDOF2][numgrids];275 276 /*Get dh1dh2dh3 in actual coordinates system : */277 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);278 279 /*Build B': */280 for (i=0;i<numgrids;i++){281 *(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i];282 *(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i];283 *(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i];284 *(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i];285 *(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i];286 *(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i];287 }288 }289 /*}}}*/290 /*FUNCTION TriaRef::GetBprimePrognostic{{{1*/291 void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss){292 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.293 * For grid i, Bi' can be expressed in the actual coordinate system294 * by:295 * Bi_prime=[ dh/dx ]296 * [ dh/dy ]297 * where h is the interpolation function for grid i.298 *299 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)300 */301 302 int i;303 const int NDOF1=1;304 const int NDOF2=2;305 const int numgrids=3;306 307 /*Same thing in the actual coordinate system: */308 double dh1dh3[NDOF2][numgrids];309 310 /*Get dh1dh2dh3 in actual coordinates system : */311 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);312 313 /*Build B': */314 for (i=0;i<numgrids;i++){315 *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i];316 *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i];317 }318 }319 /*}}}*/320 192 /*FUNCTION TriaRef::GetBprimePrognostic{{{1*/ 321 193 void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss){ … … 348 220 } 349 221 /*}}}*/ 350 /*FUNCTION TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof) {{{1*/ 351 void TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof){ 352 /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. 353 * For grid i, Li can be expressed in the actual coordinate system 354 * by: 355 * numdof=1: 356 * Li=h; 357 * numdof=2: 358 * Li=[ h 0 ] 359 * [ 0 h ] 360 * where h is the interpolation function for grid i. 361 * 362 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2) 363 */ 364 365 int i; 366 const int NDOF2=2; 367 const int numgrids=3; 368 369 double l1l2l3[3]; 370 371 372 /*Get l1l2l3 in actual coordinate system: */ 373 GetNodalFunctions(l1l2l3, gauss); 374 375 /*Build L: */ 376 if(numdof==1){ 377 for (i=0;i<numgrids;i++){ 378 L[i]=l1l2l3[i]; 379 } 380 } 381 else{ 382 for (i=0;i<numgrids;i++){ 383 *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i]; 384 *(L+numdof*numgrids*0+numdof*i+1)=0; 385 *(L+numdof*numgrids*1+numdof*i)=0; 386 *(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i]; 387 } 388 } 389 } 390 /*}}}*/ 391 /*FUNCTION TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof) {{{1*/ 222 /*FUNCTION TriaRef::GetL{{{1*/ 392 223 void TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof){ 393 224 /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. … … 428 259 } 429 260 /*}}}*/ 430 /*FUNCTION TriaRef::GetJacobian (double* J, double* xyz_list,double* gauss){{{1*/431 void TriaRef::GetJacobian(double* J, double* xyz_list, double* gauss){261 /*FUNCTION TriaRef::GetJacobian{{{1*/ 262 void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){ 432 263 /*The Jacobian is constant over the element, discard the gaussian points. 433 264 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ … … 451 282 } 452 283 /*}}}*/ 453 /*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/ 454 void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){ 455 /*The Jacobian is constant over the element, discard the gaussian points. 456 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 457 458 const int NDOF2=2; 459 const int numgrids=3; 460 double x1,y1,x2,y2,x3,y3; 461 462 x1=*(xyz_list+numgrids*0+0); 463 y1=*(xyz_list+numgrids*0+1); 464 x2=*(xyz_list+numgrids*1+0); 465 y2=*(xyz_list+numgrids*1+1); 466 x3=*(xyz_list+numgrids*2+0); 467 y3=*(xyz_list+numgrids*2+1); 468 469 470 *(J+NDOF2*0+0)=0.5*(x2-x1); 471 *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2); 472 *(J+NDOF2*0+1)=0.5*(y2-y1); 473 *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2); 474 } 475 /*}}}*/ 476 /*FUNCTION TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/ 284 /*FUNCTION TriaRef::GetSegmentJacobian{{{1*/ 477 285 void TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss){ 478 286 /*The Jacobian is constant over the element, discard the gaussian points. … … 489 297 } 490 298 /*}}}*/ 491 /*FUNCTION TriaRef::GetSegmentJacobianDeterminant (double* Jdet, double* xyz_list,GaussTria* gauss){{{1*/299 /*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{1*/ 492 300 void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){ 493 301 /*The Jacobian determinant is constant over the element, discard the gaussian points. … … 499 307 } 500 308 /*}}}*/ 501 /*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss) {{{1*/ 502 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){ 503 504 /*The Jacobian determinant is constant over the element, discard the gaussian points. 505 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 506 double J[2][2]; 507 508 /*Get Jacobian*/ 509 GetJacobian(&J[0][0],xyz_list,gauss); 510 511 /*Get Determinant*/ 512 Matrix2x2Determinant(Jdet,&J[0][0]); 513 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 514 515 } 516 /*}}}*/ 517 /*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/ 309 /*FUNCTION TriaRef::GetJacobianDeterminant2d{{{1*/ 518 310 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss){ 519 311 /*The Jacobian determinant is constant over the element, discard the gaussian points. … … 531 323 /*}}}*/ 532 324 /*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/ 533 void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list, double* gauss){325 void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss){ 534 326 /*The Jacobian determinant is constant over the element, discard the gaussian points. 535 327 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ … … 552 344 } 553 345 /*}}}*/ 554 /*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/ 555 void TriaRef::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss){ 556 /*The Jacobian determinant is constant over the element, discard the gaussian points. 557 * J is assumed to have been allocated of size NDOF2xNDOF2.*/ 558 559 double x1,x2,x3,y1,y2,y3,z1,z2,z3; 560 561 x1=*(xyz_list+3*0+0); 562 y1=*(xyz_list+3*0+1); 563 z1=*(xyz_list+3*0+2); 564 x2=*(xyz_list+3*1+0); 565 y2=*(xyz_list+3*1+1); 566 z2=*(xyz_list+3*1+2); 567 x3=*(xyz_list+3*2+0); 568 y3=*(xyz_list+3*2+1); 569 z3=*(xyz_list+3*2+2); 570 571 *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); 572 if(*Jdet<0) ISSMERROR("negative jacobian determinant!"); 573 574 } 575 /*}}}*/ 576 /*FUNCTION TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss) {{{1*/ 577 void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){ 346 /*FUNCTION TriaRef::GetJacobianInvert{{{1*/ 347 void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss){ 578 348 579 349 /*Jacobian*/ … … 588 358 } 589 359 /*}}}*/ 590 /*FUNCTION TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss) {{{1*/ 591 void TriaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss){ 592 593 /*Jacobian*/ 594 double J[2][2]; 595 596 /*Call Jacobian routine to get the jacobian:*/ 597 GetJacobian(&J[0][0], xyz_list, gauss); 598 599 /*Invert Jacobian matrix: */ 600 Matrix2x2Invert(Jinv,&J[0][0]); 601 602 } 603 /*}}}*/ 604 /*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss) {{{1*/ 605 void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){ 606 /*This routine returns the values of the nodal functions at the gaussian point.*/ 607 608 l1l2l3[0]=gauss[0]; 609 l1l2l3[1]=gauss[1]; 610 l1l2l3[2]=gauss[2]; 611 612 } 613 /*}}}*/ 614 /*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss) {{{1*/ 360 /*FUNCTION TriaRef::GetNodalFunctions{{{1*/ 615 361 void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){ 616 362 /*This routine returns the values of the nodal functions at the gaussian point.*/ … … 622 368 } 623 369 /*}}}*/ 624 /*FUNCTION TriaRef::GetSegmentNodalFunctions (double* l1l2,GaussTria* gauss){{{1*/370 /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{1*/ 625 371 void TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss,int index1,int index2){ 626 372 /*This routine returns the values of the nodal functions at the gaussian point.*/ … … 636 382 } 637 383 /*}}}*/ 638 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives (double* dh1dh3,double* xyz_list, double* gauss){{{1*/639 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){384 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{1*/ 385 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){ 640 386 641 387 /*This routine returns the values of the nodal functions derivatives (with respect to the … … 665 411 } 666 412 /*}}}*/ 667 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss) {{{1*/ 668 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){ 669 670 /*This routine returns the values of the nodal functions derivatives (with respect to the 671 * actual coordinate system): */ 672 int i; 673 const int NDOF2 = 2; 674 const int numgrids = 3; 675 double dh1dh3_ref[NDOF2][numgrids]; 676 double Jinv[NDOF2][NDOF2]; 677 678 /*Get derivative values with respect to parametric coordinate system: */ 679 GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss); 680 681 /*Get Jacobian invert: */ 682 GetJacobianInvert(&Jinv[0][0], xyz_list, gauss); 683 684 /*Build dh1dh3: 685 * 686 * [dhi/dx]= Jinv*[dhi/dr] 687 * [dhi/dy] [dhi/ds] 688 */ 689 for (i=0;i<numgrids;i++){ 690 dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i]; 691 dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i]; 692 } 693 694 } 695 /*}}}*/ 696 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss) {{{1*/ 697 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss){ 413 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{1*/ 414 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){ 698 415 699 416 /*This routine returns the values of the nodal functions derivatives (with respect to the … … 717 434 } 718 435 /*}}}*/ 719 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss) {{{1*/ 720 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){ 721 722 /*This routine returns the values of the nodal functions derivatives (with respect to the 723 * natural coordinate system) at the gaussian point. */ 724 725 const int NDOF2=2; 726 const int numgrids=3; 727 728 /*First nodal function: */ 729 *(dl1dl3+numgrids*0+0)=-0.5; 730 *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3); 731 732 /*Second nodal function: */ 733 *(dl1dl3+numgrids*0+1)=0.5; 734 *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3); 735 736 /*Third nodal function: */ 737 *(dl1dl3+numgrids*0+2)=0; 738 *(dl1dl3+numgrids*1+2)=1.0/SQRT3; 739 740 } 741 /*}}}*/ 742 /*FUNCTION TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss) {{{1*/ 743 void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss){ 436 /*FUNCTION TriaRef::GetParameterDerivativeValue{{{1*/ 437 void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss){ 744 438 745 439 /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian … … 763 457 } 764 458 /*}}}*/ 765 /*FUNCTION TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss) {{{1*/ 766 void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss){ 767 768 /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 769 * point specified by gauss_l1l2l3: 770 * dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx 771 * dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx 772 * 773 * p is a vector of size 2x1 already allocated. 774 */ 775 776 /*Nodal Derivatives*/ 777 double dh1dh3[2][3]; //nodal derivative functions in actual coordinate system. 778 779 /*Get dh1dh2dh3 in actual coordinate system: */ 780 GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss); 781 782 /*Assign values*/ 783 *(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2]; 784 *(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2]; 785 786 } 787 /*}}}*/ 788 /*FUNCTION TriaRef::GetParameterValue(double* p, double* plist, double* gauss){{{1*/ 789 void TriaRef::GetParameterValue(double* p, double* plist, double* gauss){ 459 /*FUNCTION TriaRef::GetParameterValue{{{1*/ 460 void TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){ 790 461 791 462 /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian … … 802 473 } 803 474 /*}}}*/ 804 /*FUNCTION TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){{{1*/805 void TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){806 807 /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian808 * point specifie by gauss: */809 810 /*nodal functions annd output: */811 double l1l2l3[3];812 813 /*Get nodal functions*/814 GetNodalFunctions(l1l2l3, gauss);815 816 /*Get parameter*/817 *p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2];818 }819 /*}}}*/ -
issm/trunk/src/c/objects/Elements/TriaRef.h
r5738 r5743 25 25 26 26 /*Numerics*/ 27 void GetBMacAyeal(double* B, double* xyz_list, double* gauss);28 27 void GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss); 29 void GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss);30 28 void GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss); 31 void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss);32 29 void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss); 33 void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss);34 30 void GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss); 35 void GetL(double* L, double* xyz_list,double* gauss,int numdof);36 31 void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof); 37 void GetJacobian(double* J, double* xyz_list,double* gauss);38 32 void GetJacobian(double* J, double* xyz_list,GaussTria* gauss); 39 33 void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss); 40 34 void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss); 41 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);42 35 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss); 43 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);44 36 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss); 45 void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss);46 37 void GetJacobianInvert(double* Jinv, double* xyz_list,GaussTria* gauss); 47 void GetNodalFunctions(double* l1l2l3, double* gauss);48 38 void GetNodalFunctions(double* l1l2l3,GaussTria* gauss); 49 39 void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2); 50 40 void GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2); 51 41 void GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2); 52 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);53 42 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss); 54 void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss);55 43 void GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss); 56 void GetParameterValue(double* pp, double* plist, double* gauss);57 44 void GetParameterValue(double* pp, double* plist, GaussTria* gauss); 58 void GetParameterDerivativeValue(double* pp, double* plist,double* xyz_list, double* gauss);59 45 void GetParameterDerivativeValue(double* pp, double* plist,double* xyz_list, GaussTria* gauss); 60 46 -
issm/trunk/src/c/objects/Inputs/BoolInput.cpp
r5660 r5743 166 166 void BoolInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");} 167 167 /*}}}*/ 168 /*FUNCTION BoolInput::GetParameterValue(double* pvalue,double* gauss){{{1*/169 void BoolInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}170 /*}}}*/171 168 /*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ 172 169 void BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} … … 174 171 /*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 175 172 void BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} 176 /*}}}*/177 /*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/178 void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}179 173 /*}}}*/ 180 174 /*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ -
issm/trunk/src/c/objects/Inputs/BoolInput.h
r5736 r5743 47 47 void GetParameterValue(int* pvalue); 48 48 void GetParameterValue(double* pvalue); 49 void GetParameterValue(double* pvalue,double* gauss);50 49 void GetParameterValue(double* pvalue,GaussTria* gauss); 51 50 void GetParameterValue(double* pvalue,GaussPenta* gauss); 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);53 51 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 54 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); -
issm/trunk/src/c/objects/Inputs/DoubleInput.cpp
r5707 r5743 182 182 } 183 183 /*}}}*/ 184 /*FUNCTION DoubleInput::GetParameterValue(double* pvalue,double* gauss){{{1*/185 void DoubleInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}186 /*}}}*/187 184 /*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ 188 185 void DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} … … 190 187 /*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 191 188 void DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} 192 /*}}}*/193 /*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/194 void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}195 189 /*}}}*/ 196 190 /*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ -
issm/trunk/src/c/objects/Inputs/DoubleInput.h
r5736 r5743 46 46 void GetParameterValue(int* pvalue); 47 47 void GetParameterValue(double* pvalue); 48 void GetParameterValue(double* pvalue,double* gauss);49 48 void GetParameterValue(double* pvalue,GaussTria* gauss); 50 49 void GetParameterValue(double* pvalue,GaussPenta* gauss); 51 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);52 50 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 53 51 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); -
issm/trunk/src/c/objects/Inputs/Input.h
r5736 r5743 25 25 virtual void GetParameterValue(int* pvalue)=0; 26 26 virtual void GetParameterValue(double* pvalue)=0; 27 virtual void GetParameterValue(double* pvalue,double* gauss)=0;28 27 virtual void GetParameterValue(double* pvalue,GaussTria* gauss)=0; 29 28 virtual void GetParameterValue(double* pvalue,GaussPenta* gauss)=0; 30 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss)=0;31 29 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss)=0; 32 30 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss)=0; -
issm/trunk/src/c/objects/Inputs/IntInput.cpp
r5660 r5743 167 167 } 168 168 /*}}}*/ 169 /*FUNCTION IntInput::GetParameterValue(double* pvalue,double* gauss){{{1*/170 void IntInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}171 /*}}}*/172 169 /*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ 173 170 void IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");} … … 175 172 /*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 176 173 void IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");} 177 /*}}}*/178 /*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/179 void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}180 174 /*}}}*/ 181 175 /*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/ -
issm/trunk/src/c/objects/Inputs/IntInput.h
r5736 r5743 47 47 void GetParameterValue(int* pvalue); 48 48 void GetParameterValue(double* pvalue); 49 void GetParameterValue(double* pvalue,double* gauss);50 49 void GetParameterValue(double* pvalue,GaussTria* gauss); 51 50 void GetParameterValue(double* pvalue,GaussPenta* gauss); 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);53 51 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 54 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); -
issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp
r5736 r5743 176 176 177 177 /*Object functions*/ 178 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue, double* gauss){{{1*/179 void PentaVertexInput::GetParameterValue(double* pvalue, double* gauss){178 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/ 179 void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ 180 180 181 181 /*Call PentaRef function*/ 182 182 PentaRef::GetParameterValue(pvalue,&values[0],gauss); 183 183 184 }185 /*}}}*/186 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/187 void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){188 189 /*Call PentaRef function*/190 PentaRef::GetParameterValue(pvalue,&values[0],gauss);191 192 }193 /*}}}*/194 /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/195 void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){196 197 /*Call PentaRef function*/198 PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);199 184 } 200 185 /*}}}*/ -
issm/trunk/src/c/objects/Inputs/PentaVertexInput.h
r5736 r5743 47 47 void GetParameterValue(int* pvalue){ISSMERROR("not implemented yet");}; 48 48 void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");}; 49 void GetParameterValue(double* pvalue,double* gauss);50 49 void GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR("not implemented yet");}; 51 50 void GetParameterValue(double* pvalue,GaussPenta* gauss); 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);53 51 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");}; 54 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss); -
issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
r5736 r5743 165 165 166 166 /*Object functions*/ 167 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue, double* gauss){{{1*/168 void TriaVertexInput::GetParameterValue(double* pvalue, double* gauss){167 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/ 168 void TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){ 169 169 170 170 /*Call TriaRef function*/ 171 171 TriaRef::GetParameterValue(pvalue,&values[0],gauss); 172 172 173 }174 /*}}}*/175 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/176 void TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){177 178 /*Call TriaRef function*/179 TriaRef::GetParameterValue(pvalue,&values[0],gauss);180 181 }182 /*}}}*/183 /*FUNCTION TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/184 void TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){185 186 /*Call TriaRef function*/187 TriaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);188 173 } 189 174 /*}}}*/ -
issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
r5736 r5743 47 47 void GetParameterValue(int* pvalue){ISSMERROR("not implemented yet");} 48 48 void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");} 49 void GetParameterValue(double* pvalue,double* gauss);50 49 void GetParameterValue(double* pvalue,GaussTria* gauss); 51 50 void GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR("not implemented yet");}; 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);53 51 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss); 54 52 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");}; -
issm/trunk/src/c/objects/Loads/Friction.cpp
r5682 r5743 55 55 } 56 56 /*}}}*/ 57 /*FUNCTION Friction::GetAlpha2 {{{1*/58 void Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){57 /*FUNCTION Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{1*/ 58 void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){ 59 59 60 60 /*This routine calculates the basal friction coefficient 61 alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/61 alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ 62 62 63 63 /*diverse: */ … … 87 87 r=drag_q/drag_p; 88 88 s=1./drag_p; 89 89 90 90 //From bed and thickness, compute effective pressure when drag is viscous: 91 91 Neff=gravity*(rho_ice*thickness+rho_water*bed); 92 92 93 93 /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 94 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour95 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should96 replace it by Neff=0 (ie, equival it to an ice shelf)*/94 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 95 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 96 replace it by Neff=0 (ie, equival it to an ice shelf)*/ 97 97 if (Neff<0)Neff=0; 98 98 … … 109 109 } 110 110 else ISSMERROR("element_type %s not supported yet",element_type); 111 111 112 112 alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); 113 113 … … 116 116 } 117 117 /*}}}*/ 118 /*FUNCTION Friction::GetAlpha2 {{{1*/119 void Friction::GetAlpha2(double* palpha2, Gauss Tria* gauss,int vxenum,int vyenum,int vzenum){118 /*FUNCTION Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{1*/ 119 void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){ 120 120 121 121 /*This routine calculates the basal friction coefficient … … 177 177 } 178 178 /*}}}*/ 179 /*FUNCTION Friction::GetAlpha2{{{1*/180 void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){181 182 /*This routine calculates the basal friction coefficient183 alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/184 185 /*diverse: */186 double r,s;187 double drag_p, drag_q;188 double gravity,rho_ice,rho_water;189 double Neff;190 double thickness,bed;191 double vx,vy,vz,vmag;192 double drag_coefficient;193 double alpha2;194 195 196 /*Recover parameters: */197 inputs->GetParameterValue(&drag_p,DragPEnum);198 inputs->GetParameterValue(&drag_q,DragQEnum);199 inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);200 inputs->GetParameterValue(&bed, gauss,BedEnum);201 inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);202 203 /*Get material parameters: */204 gravity=matpar->GetG();205 rho_ice=matpar->GetRhoIce();206 rho_water=matpar->GetRhoWater();207 208 //compute r and q coefficients: */209 r=drag_q/drag_p;210 s=1./drag_p;211 212 //From bed and thickness, compute effective pressure when drag is viscous:213 Neff=gravity*(rho_ice*thickness+rho_water*bed);214 215 /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because216 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour217 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should218 replace it by Neff=0 (ie, equival it to an ice shelf)*/219 if (Neff<0)Neff=0;220 221 if(strcmp(element_type,"2d")==0){222 inputs->GetParameterValue(&vx, gauss,vxenum);223 inputs->GetParameterValue(&vy, gauss,vyenum);224 vmag=sqrt(pow(vx,2)+pow(vy,2));225 }226 else if (strcmp(element_type,"3d")==0){227 inputs->GetParameterValue(&vx, gauss,vxenum);228 inputs->GetParameterValue(&vy, gauss,vyenum);229 inputs->GetParameterValue(&vz, gauss,vzenum);230 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));231 }232 else ISSMERROR("element_type %s not supported yet",element_type);233 234 alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));235 236 /*Assign output pointers:*/237 *palpha2=alpha2;238 }239 /*}}}*/240 179 /*FUNCTION Friction::GetAlphaComplement {{{1*/ 241 void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){242 180 void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){ 181 243 182 /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 244 183 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: … … 274 213 r=drag_q/drag_p; 275 214 s=1./drag_p; 276 215 277 216 //From bed and thickness, compute effective pressure when drag is viscous: 278 217 Neff=gravity*(rho_ice*thickness+rho_water*bed); … … 288 227 inputs->GetParameterValue(&vy, gauss,vyenum); 289 228 vmag=sqrt(pow(vx,2)+pow(vy,2)); 290 229 291 230 alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); 292 231 … … 296 235 } 297 236 /*}}}*/ 298 /*FUNCTION Friction::GetAlphaComplement {{{1*/299 void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){300 301 /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.302 * FrictionGetAlphaComplement is used in control methods on drag, and it computes:303 * alpha_complement= Neff ^r * vel ^s*/304 305 /*diverse: */306 int i;307 double Neff;308 double r,s;309 double vx;310 double vy;311 double vmag;312 double drag_p,drag_q;313 double drag_coefficient;314 double bed,thickness;315 double gravity,rho_ice,rho_water;316 double alpha_complement;317 318 /*Recover parameters: */319 inputs->GetParameterValue(&drag_p,DragPEnum);320 inputs->GetParameterValue(&drag_q,DragQEnum);321 inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);322 inputs->GetParameterValue(&bed, gauss,BedEnum);323 inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);324 325 /*Get material parameters: */326 gravity=matpar->GetG();327 rho_ice=matpar->GetRhoIce();328 rho_water=matpar->GetRhoWater();329 330 331 //compute r and q coefficients: */332 r=drag_q/drag_p;333 s=1./drag_p;334 335 //From bed and thickness, compute effective pressure when drag is viscous:336 Neff=gravity*(rho_ice*thickness+rho_water*bed);337 338 /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because339 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour340 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should341 replace it by Neff=0 (ie, equival it to an ice shelf)*/342 if (Neff<0)Neff=0;343 344 //We need the velocity magnitude to evaluate the basal stress:345 inputs->GetParameterValue(&vx, gauss,vxenum);346 inputs->GetParameterValue(&vy, gauss,vyenum);347 vmag=sqrt(pow(vx,2)+pow(vy,2));348 349 alpha_complement=pow(Neff,r)*pow(vmag,(s-1));350 351 /*Assign output pointers:*/352 *palpha_complement=alpha_complement;353 354 }355 /*}}}*/ -
issm/trunk/src/c/objects/Loads/Friction.h
r5682 r5743 27 27 28 28 void Echo(void); 29 void GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);30 29 void GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum); 31 30 void GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); 32 void GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);33 31 void GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum); 34 32 -
issm/trunk/src/c/objects/Loads/Icefront.cpp
r5742 r5743 549 549 /*Objects: */ 550 550 double pe_g[numdofs]={0.0}; 551 Node** element_nodes=NULL;552 551 int element_type; 553 552 … … 579 578 580 579 //Identify which grids are comprised in the quad: 581 if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*)); 582 element->GetNodes((void**)element_nodes); 583 584 grid1=-1; grid2=-1; grid3=-1; grid4=-1; 585 for(i=0;i<6;i++){ 586 if (nodes[0]==element_nodes[i])grid1=i; 587 if (nodes[1]==element_nodes[i])grid2=i; 588 if (nodes[2]==element_nodes[i])grid3=i; 589 if (nodes[3]==element_nodes[i])grid4=i; 590 } 591 592 if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){ 593 ISSMERROR("could not find element grids corresponding to quad icefront!"); 594 } 580 grid1=element->GetNodeIndex(nodes[0]); 581 grid2=element->GetNodeIndex(nodes[1]); 582 grid3=element->GetNodeIndex(nodes[2]); 583 grid4=element->GetNodeIndex(nodes[3]); 595 584 596 585 /*Build new xyz, bed and thickness lists for grids 1 to 4: */ … … 653 642 654 643 /*Free ressources:*/ 655 xfree((void**)&element_nodes);656 644 xfree((void**)&doflist); 657 645 … … 678 666 /*Objects: */ 679 667 double pe_g[numdofs]={0.0}; 680 Node** element_nodes=NULL;681 668 Penta* penta=NULL; 682 669 … … 709 696 710 697 //Identify which grids are comprised in the quad: 711 element_nodes=(Node**)xmalloc(6*sizeof(Node*)); 712 element->GetNodes((void**)element_nodes); 713 714 grid1=-1; grid2=-1; grid3=-1; grid4=-1; 715 for(i=0;i<6;i++){ 716 if (nodes[0]==element_nodes[i])grid1=i; 717 if (nodes[1]==element_nodes[i])grid2=i; 718 if (nodes[2]==element_nodes[i])grid3=i; 719 if (nodes[3]==element_nodes[i])grid4=i; 720 } 721 722 if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){ 723 ISSMERROR("could not find element grids corresponding to quad icefront!"); 724 } 698 grid1=element->GetNodeIndex(nodes[0]); 699 grid2=element->GetNodeIndex(nodes[1]); 700 grid3=element->GetNodeIndex(nodes[2]); 701 grid4=element->GetNodeIndex(nodes[3]); 725 702 726 703 /*Build new xyz, bed and thickness lists for grids 1 to 4: */ … … 784 761 785 762 /*Free ressources:*/ 786 xfree((void**)&element_nodes);787 763 xfree((void**)&doflist); 788 764 … … 1019 995 complete_list[j][2]=0; 1020 996 } 1021 tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0], l1l2l3_tria);997 tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss); 1022 998 } 1023 999 … … 1261 1237 complete_list[j][2]=0; 1262 1238 } 1263 tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0], l1l2l3_tria);1239 tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss); 1264 1240 } 1265 1241 … … 1310 1286 1311 1287 } //for(i=0;i<4;i++) 1312 } //for(ig=0;ig<num_gauss;ig++)1288 } 1313 1289 1314 1290 delete tria;
Note:
See TracChangeset
for help on using the changeset viewer.