Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6722) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6723) @@ -699,4 +699,39 @@ } /*}}}*/ +/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokes{{{1*/ +ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokes(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); + ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + _error_("KMatrix coupling MacAyealStokes not terminated yet"); + ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2); + delete Ke1; + delete Ke2; + Ke1=CreateKMatrixDiagnosticMacAyeal2d(); + Ke2=CreateKMatrixDiagnosticStokes(); + + + /*Constants*/ + const int numdofm=NDOF2*NUMVERTICES; + const int numdofs=NDOF4*NUMVERTICES; + const int numdoftotal=(NDOF2+NDOF4)*NUMVERTICES; + int i,j; + + for(i=0;ivalues[(i+numdofm)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF4*j+0]; + Ke->values[(i+numdofm)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF4*j+1]; + } + for(i=0;ivalues[i*numdoftotal+numdofm+NDOF4*j+0]+=Ke1->values[i*numdofm+NDOF2*j+0]; + Ke->values[i*numdoftotal+numdofm+NDOF4*j+1]+=Ke1->values[i*numdofm+NDOF2*j+1]; + } + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ /*FUNCTION Penta::CreateKMatrixCouplingPattynStokes{{{1*/ ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){ @@ -751,4 +786,6 @@ case MacAyealPattynApproximationEnum: return CreateKMatrixDiagnosticMacAyealPattyn(); + case MacAyealStokesApproximationEnum: + return CreateKMatrixDiagnosticMacAyealStokes(); case PattynStokesApproximationEnum: return CreateKMatrixDiagnosticPattynStokes(); @@ -967,4 +1004,20 @@ } /*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealStokes{{{1*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealStokes(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal2d(); + ElementMatrix* Ke2=CreateKMatrixDiagnosticStokes(); + ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealStokes(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + delete Ke3; + return Ke; +} +/*}}}*/ /*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{1*/ ElementMatrix* Penta::CreateKMatrixDiagnosticPattyn(void){ @@ -1112,5 +1165,5 @@ /*If on water or not Stokes, skip stiffness: */ inputs->GetParameterValue(&approximation,ApproximationEnum); - if(IsOnWater() || (approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; + if(IsOnWater() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); @@ -1180,5 +1233,5 @@ /*If on water or not Stokes, skip stiffness: */ inputs->GetParameterValue(&approximation,ApproximationEnum); - if(IsOnWater() || IsOnShelf() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; + if(IsOnWater() || IsOnShelf() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum)) return NULL; ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); @@ -1696,10 +1749,11 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/ -ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){ - +/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokes {{{1*/ +ElementVector* Penta::CreatePVectorCouplingMacAyealStokes(void){ + + _error_("coupling MacAyeal Stokes not implemented yet"); /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous(); - ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction(); + ElementVector* pe1=CreatePVectorCouplingMacAyealStokesViscous(); + ElementVector* pe2=CreatePVectorCouplingMacAyealStokesFriction(); ElementVector* pe =new ElementVector(pe1,pe2); @@ -1710,6 +1764,6 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{1*/ -ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){ +/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesViscous {{{1*/ +ElementVector* Penta::CreatePVectorCouplingMacAyealStokesViscous(void){ /*Constants*/ @@ -1731,5 +1785,5 @@ if(IsOnWater()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=PattynStokesApproximationEnum) return NULL; + if(approximation!=MacAyealStokesApproximationEnum) return NULL; ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); @@ -1740,5 +1794,5 @@ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); _assert_(vzpattyn_input); + Input* vzpattyn_input=inputs->GetInput(VzMacAyealEnum); _assert_(vzpattyn_input); /* Start looping on the number of gaussian points: */ @@ -1770,6 +1824,6 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/ -ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){ +/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesFriction{{{1*/ +ElementVector* Penta::CreatePVectorCouplingMacAyealStokesFriction(void){ /*Constants*/ @@ -1795,4 +1849,149 @@ if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation!=MacAyealStokesApproximationEnum) return NULL; + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* vzpattyn_input=inputs->GetInput(VzMacAyealEnum); _assert_(vzpattyn_input); + + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(l1l6, gauss); + + vzpattyn_input->GetParameterValue(&w, gauss); + vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); + + BedNormal(&bed_normal[0],xyz_list_tria); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + + for(i=0;ivalues[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*l1l6[i]; + pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*l1l6[i]; + pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i]; + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/ +ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous(); + ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction(); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{1*/ +ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + + /*Intermediaries */ + int i,j,ig; + int approximation; + double viscosity,Jdet; + double stokesreconditioning; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double dw[3]; + double xyz_list[NUMVERTICES][3]; + double l1l6[6]; //for the six nodes of the penta + double dh1dh6[3][6]; //for the six nodes of the penta + GaussPenta *gauss=NULL; + + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation!=PattynStokesApproximationEnum) return NULL; + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); _assert_(vzpattyn_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); + + vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + for(i=0;ivalues[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]; + pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]; + pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+2*dw[2]*dh1dh6[2][i]); + pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i]; + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/ +ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF4; + + /*Intermediaries*/ + int i,j,ig; + int approximation,analysis_type; + double Jdet,Jdet2d; + double stokesreconditioning; + double bed_normal[3]; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity, w, alpha2_gauss; + double dw[3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double xyz_list[NUMVERTICES][3]; + double l1l6[6]; //for the six nodes of the penta + Tria* tria=NULL; + Friction* friction=NULL; + GaussPenta *gauss=NULL; + + /*Initialize Element vector and return if necessary*/ + if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL; + inputs->GetParameterValue(&approximation,ApproximationEnum); if(approximation!=PattynStokesApproximationEnum) return NULL; ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); @@ -1860,4 +2059,6 @@ case MacAyealPattynApproximationEnum: return CreatePVectorDiagnosticMacAyealPattyn(); + case MacAyealStokesApproximationEnum: + return CreatePVectorDiagnosticMacAyealStokes(); case PattynStokesApproximationEnum: return CreatePVectorDiagnosticPattynStokes(); @@ -1878,4 +2079,20 @@ delete pe1; delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealStokes{{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticMacAyealStokes(void){ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorDiagnosticMacAyeal(); + ElementVector* pe2=CreatePVectorDiagnosticStokes(); + ElementVector* pe3=CreatePVectorCouplingMacAyealStokes(); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; return pe; } @@ -2087,5 +2304,5 @@ if(IsOnWater()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; + if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); @@ -2157,5 +2374,5 @@ if(IsOnWater() || !IsOnBed() || !IsOnShelf()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); - if(approximation!=StokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; + if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL; ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); @@ -3603,4 +3820,7 @@ else if (*(iomodel->elements_type+index)==StokesApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==MacAyealStokesApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum)); } else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 6722) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 6723) @@ -131,4 +131,5 @@ ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void); ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void); + ElementMatrix* CreateKMatrixCouplingMacAyealStokes(void); ElementMatrix* CreateKMatrixCouplingPattynStokes(void); ElementMatrix* CreateKMatrixDiagnosticHoriz(void); @@ -139,4 +140,5 @@ ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void); ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void); + ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void); ElementMatrix* CreateKMatrixDiagnosticPattyn(void); ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void); @@ -160,4 +162,7 @@ ElementVector* CreatePVectorAdjointMacAyeal(void); ElementVector* CreatePVectorAdjointPattyn(void); + ElementVector* CreatePVectorCouplingMacAyealStokes(void); + ElementVector* CreatePVectorCouplingMacAyealStokesViscous(void); + ElementVector* CreatePVectorCouplingMacAyealStokesFriction(void); ElementVector* CreatePVectorCouplingPattynStokes(void); ElementVector* CreatePVectorCouplingPattynStokesViscous(void); @@ -167,4 +172,5 @@ ElementVector* CreatePVectorDiagnosticMacAyeal(void); ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void); + ElementVector* CreatePVectorDiagnosticMacAyealStokes(void); ElementVector* CreatePVectorDiagnosticPattyn(void); ElementVector* CreatePVectorDiagnosticPattynStokes(void);