Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5925) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5926) @@ -740,4 +740,5 @@ /*retrive parameters: */ + ElementVector* pe=NULL; int analysis_type; parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -750,8 +751,12 @@ switch(analysis_type){ case DiagnosticHorizAnalysisEnum: - CreatePVectorDiagnosticHoriz( pg); + pe=CreatePVectorDiagnosticHoriz(); + if(pe) pe->AddToGlobal(pg,pf); + delete pe; break; case AdjointHorizAnalysisEnum: - CreatePVectorAdjointHoriz( pg); + pe=CreatePVectorAdjointHoriz(); + if(pe) pe->AddToGlobal(pg,pf); + delete pe; break; case DiagnosticHutterAnalysisEnum: @@ -2945,5 +2950,5 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/ -void Penta::CreatePVectorAdjointHoriz(Vec pg){ +ElementVector* Penta::CreatePVectorAdjointHoriz(void){ int approximation; @@ -2952,14 +2957,11 @@ switch(approximation){ case MacAyealApproximationEnum: - CreatePVectorAdjointMacAyeal( pg); - break; + return CreatePVectorAdjointMacAyeal(); case PattynApproximationEnum: - CreatePVectorAdjointPattyn( pg); - break; + return CreatePVectorAdjointPattyn(); case NoneApproximationEnum: - return; + return NULL; case StokesApproximationEnum: - CreatePVectorAdjointStokes( pg); - break; + return CreatePVectorAdjointStokes(); default: ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); @@ -2968,7 +2970,7 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorAdjointMacAyeal{{{1*/ -void Penta::CreatePVectorAdjointMacAyeal(Vec p_g){ - - if (!IsOnBed() || IsOnWater()) return; +ElementVector* Penta::CreatePVectorAdjointMacAyeal(){ + + if (!IsOnBed() || IsOnWater()) return NULL; /*Call Tria function*/ @@ -2976,15 +2978,13 @@ ElementVector* pe=tria->CreatePVectorAdjointHoriz(); delete tria->matice; delete tria; - pe->AddToGlobal(p_g,NULL); - delete pe; /*clean up and return*/ - return; + return pe; } /*}}}*/ /*FUNCTION Penta::CreatePVectorAdjointPattyn{{{1*/ -void Penta::CreatePVectorAdjointPattyn(Vec p_g){ - - if (!IsOnSurface() || IsOnWater()) return; +ElementVector* Penta::CreatePVectorAdjointPattyn(void){ + + if (!IsOnSurface() || IsOnWater()) return NULL; /*Call Tria function*/ @@ -2992,9 +2992,7 @@ ElementVector* pe=tria->CreatePVectorAdjointHoriz(); delete tria->matice; delete tria; - pe->AddToGlobal(p_g,NULL); - delete pe; /*clean up and return*/ - return; + return pe; } /*}}}*/ @@ -3048,25 +3046,30 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{1*/ -void Penta::CreatePVectorCouplingPattynStokes( Vec pg){ - - /*indexing: */ +ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){ + + /*compute all load vectorsfor 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){ + + const int numdof=NUMVERTICES*NDOF4; int i,j; - const int numdof=NUMVERTICES*NDOF4; - int* doflist=NULL; - - /*parameters: */ + int ig; double xyz_list_tria[NUMVERTICES2D][3]; double xyz_list[NUMVERTICES][3]; double bed_normal[3]; - - /* gaussian points: */ - int ig; GaussPenta *gauss=NULL; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ double viscosity, w, alpha2_gauss; double dw[3]; - - /*matrices: */ double Pe_gaussian[24]={0.0}; //for the six nodes double l1l6[6]; //for the six nodes of the penta @@ -3074,28 +3077,20 @@ double Jdet; double Jdet2d; - Tria* tria=NULL; Friction* friction=NULL; - - /*parameters: */ double stokesreconditioning; int approximation; int analysis_type; - /*retrieve inputs :*/ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*retrieve some parameters: */ + if(approximation!=PattynStokesApproximationEnum) return NULL; + ElementVector* pe=NewElementVector(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); - - /*If on water or not Stokes, skip load: */ - if(IsOnWater() || approximation!=PattynStokesApproximationEnum) return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,StokesApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -3109,86 +3104,116 @@ gauss->GaussPoint(ig); - /*Compute strain rate and viscosity: */ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - /* Get nodal functions */ GetNodalFunctionsP1(&l1l6[0], gauss); GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); - /*Compute the derivative of VzPattyn*/ vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - /* Build gaussian vector */ for(i=0;iweight*viscosity*dw[0]*dh1dh6[2][i]/2; - Pe_gaussian[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2; - Pe_gaussian[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2; - Pe_gaussian[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i]; - } - } + pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2; + pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2; + pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2; + pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i]; + } + } + + /*Clean up and return*/ delete gauss; - - /*Deal with 2d friction at the bedrock interface: */ - if(IsOnBed() && !IsOnShelf()){ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{1*/ +ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){ + + const int numdof=NUMVERTICES*NDOF4; + int i,j; + int ig; + double xyz_list_tria[NUMVERTICES2D][3]; + double xyz_list[NUMVERTICES][3]; + double bed_normal[3]; + GaussPenta *gauss=NULL; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity, w, alpha2_gauss; + double dw[3]; + double Pe_gaussian[24]={0.0}; //for the six nodes + double l1l6[6]; //for the six nodes of the penta + double dh1dh6[3][6]; //for the six nodes of the penta + double Jdet; + double Jdet2d; + Tria* tria=NULL; + Friction* friction=NULL; + double stokesreconditioning; + int approximation; + int analysis_type; + + /*Initialize Element vector and return if necessary*/ + if(IsOnWater() || !IsOnBed() || IsOnShelf()) return NULL; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation!=PattynStokesApproximationEnum) return NULL; + ElementVector* pe=NewElementVector(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); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); + Input* vzpattyn_input=inputs->GetInput(VzPattynEnum); ISSMASSERT(vzpattyn_input); + + + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Get the Jacobian determinant */ + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + + /*Get L matrix : */ + GetNodalFunctionsP1(l1l6, gauss); + + /*Get normal vecyor to the bed */ + BedNormal(&bed_normal[0],xyz_list_tria); + + /*Get Viscosity*/ + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /*Get friction*/ + friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + + /*Get w and its derivatives*/ + vzpattyn_input->GetParameterValue(&w, gauss); + vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Get the Jacobian determinant */ - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - - /*Get L matrix : */ - GetNodalFunctionsP1(l1l6, gauss); - - /*Get normal vecyor to the bed */ - BedNormal(&bed_normal[0],xyz_list_tria); - - /*Get Viscosity*/ - this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /*Get friction*/ - friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); - - /*Get w and its derivatives*/ - vzpattyn_input->GetParameterValue(&w, gauss); - vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss); - - for(i=0;iweight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+viscosity*dw[2]*bed_normal[0])*l1l6[i]; - Pe_gaussian[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+viscosity*dw[2]*bed_normal[1])*l1l6[i]; - Pe_gaussian[i*NDOF4+2]+=Jdet2d*gauss->weight*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i]; - } - } - /*Free ressources:*/ - delete gauss; - } //if ( (IsOnBed()==1) && (IsOnShelf()==1)) - - /*Add Pe_reduced to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)Pe_gaussian,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&doflist); - + pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[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]+viscosity*dw[2]*bed_normal[1])*l1l6[i]; + pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*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::CreatePVectorDiagnosticHoriz{{{1*/ -void Penta::CreatePVectorDiagnosticHoriz(Vec pg){ +ElementVector* Penta::CreatePVectorDiagnosticHoriz(void){ int approximation; @@ -3197,28 +3222,50 @@ switch(approximation){ case MacAyealApproximationEnum: - CreatePVectorDiagnosticMacAyeal( pg); - break; + return CreatePVectorDiagnosticMacAyeal(); case PattynApproximationEnum: - CreatePVectorDiagnosticPattyn( pg); - break; + return CreatePVectorDiagnosticPattyn(); case HutterApproximationEnum: - return; + return NULL; case NoneApproximationEnum: - return; + return NULL; case StokesApproximationEnum: - CreatePVectorDiagnosticStokes( pg); - break; + return CreatePVectorDiagnosticStokes(); case MacAyealPattynApproximationEnum: - CreatePVectorDiagnosticMacAyeal( pg); - CreatePVectorDiagnosticPattyn( pg); - break; + return CreatePVectorDiagnosticMacAyealPattyn(); case PattynStokesApproximationEnum: - CreatePVectorDiagnosticPattyn( pg); - CreatePVectorDiagnosticStokes( pg); - CreatePVectorCouplingPattynStokes( pg); - break; + return CreatePVectorDiagnosticPattynStokes(); default: ISSMERROR("Approximation %s not supported yet",EnumToString(approximation)); } +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealPattyn{{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticMacAyealPattyn(void){ + + /*compute all load vectorsfor this element*/ + ElementVector* pe1=CreatePVectorDiagnosticMacAyeal(); + ElementVector* pe2=CreatePVectorDiagnosticPattyn(); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticPattynStokes{{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticPattynStokes(void){ + + /*compute all load vectorsfor this element*/ + ElementVector* pe1=CreatePVectorDiagnosticPattyn(); + ElementVector* pe2=CreatePVectorDiagnosticStokes(); + ElementVector* pe3=CreatePVectorCouplingPattynStokes(); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; } /*}}}*/ @@ -3329,10 +3376,7 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorDiagnosticMacAyeal{{{1*/ -void Penta::CreatePVectorDiagnosticMacAyeal(Vec pg){ - - /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the - bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build - the stiffness matrix. */ - if (!IsOnBed() || IsOnWater()) return; +ElementVector* Penta::CreatePVectorDiagnosticMacAyeal(void){ + + if (!IsOnBed() || IsOnWater()) return NULL; /*Call Tria function*/ @@ -3340,50 +3384,31 @@ ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); delete tria->matice; delete tria; - pe->AddToGlobal(pg,NULL); - delete pe; - - /*clean up and return*/ - return; + + /*Clean up and return*/ + return pe; } /*}}}*/ /*FUNCTION Penta::CreatePVectorDiagnosticPattyn{{{1*/ -void Penta::CreatePVectorDiagnosticPattyn( Vec pg){ - - int i,j; +ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){ /* node data: */ const int numdof=NDOF2*NUMVERTICES; + int i,j; + int ig; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* parameters: */ double slope[3]; //do not put 2! this goes into GetParameterDerivativeValue, which addresses slope[3] also! double driving_stress_baseline; double thickness; - - /* gaussian points: */ - int ig; GaussPenta *gauss=NULL; - - /* Jacobian: */ double Jdet; - - /*nodal functions: */ double l1l6[6]; - - /*element vector at the gaussian points: */ - double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; - /*If on water, skip load: */ - if(IsOnWater())return; - - /*Implement standard penta element: */ - - /* Get node coordinates and dof list: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,PattynApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); @@ -3395,64 +3420,50 @@ gauss->GaussPoint(ig); - /*Compute thickness at gaussian point: */ thickness_input->GetParameterValue(&thickness, gauss); - - /*Compute slope at gaussian point: */ surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); - /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get nodal functions: */ GetNodalFunctionsP1(l1l6, gauss); - /*Compute driving stress: */ driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); - /*Build pe_g_gaussian vector: */ - for (i=0;iweight*l1l6[i]; - } - } - - /*Add pe_g_gaussian vector to pe_g: */ - for( i=0; ivalues[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i]; + } + + /*Clean up and return*/ delete gauss; + return pe; } /*}}}*/ /*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{1*/ -void Penta::CreatePVectorDiagnosticStokes( Vec pg){ - - /*indexing: */ +ElementVector* Penta::CreatePVectorDiagnosticStokes(void){ + + /*compute all load vectorsfor this element*/ + ElementVector* pe1=CreatePVectorDiagnosticStokesViscous(); + ElementVector* pe2=CreatePVectorDiagnosticStokesShelf(); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticStokesViscous {{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticStokesViscous(void){ + int i,j; + int ig; const int numdof=NUMVERTICES*NDOF4; int numdof2d=NUMVERTICES2D*NDOF4; - int* doflist=NULL; - - /*Material properties: */ double gravity,rho_ice,rho_water; - - /*parameters: */ double xyz_list_tria[NUMVERTICES2D][3]; double xyz_list[NUMVERTICES][3]; double bed_normal[3]; double bed; - - /* gaussian points: */ - int ig; GaussPenta *gauss=NULL; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ double viscosity; double water_pressure; - - /*matrices: */ double Pe_temp[27]={0.0}; //for the six nodes and the bubble double Pe_gaussian[27]={0.0}; //for the six nodes and the bubble @@ -3470,32 +3481,20 @@ double D_scalar; double tBD[27][8]; - Tria* tria=NULL; - - /*parameters: */ double stokesreconditioning; - - /*inputs: */ int approximation; - /*retrieve inputs :*/ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*retrieve some parameters: */ + if(approximation!=StokesApproximationEnum) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - - /*If on water or not Stokes, skip load: */ - if(IsOnWater() || approximation!=StokesApproximationEnum) return; - - /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); - - /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,StokesApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -3509,12 +3508,7 @@ gauss->GaussPoint(ig); - /*Compute strain rate and viscosity: */ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions */ GetNodalFunctionsMINI(&l1l7[0], gauss); @@ -3547,70 +3541,106 @@ for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; } + + /*Condensation*/ + ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_temp[0]); + + /*Clean up and return*/ delete gauss; - - /*Deal with 2d friction at the bedrock interface: */ - if(IsOnBed() && IsOnShelf()){ - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - /*Get the Jacobian determinant */ - GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); - - /* Get bed at gaussian point */ - bed_input->GetParameterValue(&bed, gauss); - - /*Get L matrix : */ - GetNodalFunctionsP1(l1l6, gauss); - - /*Get water_pressure at gaussian point */ - water_pressure=gravity*rho_water*bed; - - /*Get normal vecyor to the bed */ - BedNormal(&bed_normal[0],xyz_list_tria); - - for(i=0;iweight*Jdet2d*l1l6[i]*bed_normal[j]; - } - } - } - /*Free ressources:*/ - delete gauss; - } //if ( (IsOnBed()==1) && (IsOnShelf()==1)) - - /*Reduce the matrix */ - ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]); - - /*Add Pe_reduced to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)Pe_reduced,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&doflist); - + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticStokesShelf{{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticStokesShelf(void){ + + int i,j; + int ig; + const int numdof=NUMVERTICES*NDOF4; + int numdof2d=NUMVERTICES2D*NDOF4; + double gravity,rho_ice,rho_water; + double xyz_list_tria[NUMVERTICES2D][3]; + double xyz_list[NUMVERTICES][3]; + double bed_normal[3]; + double bed; + GaussPenta *gauss=NULL; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity; + double water_pressure; + double Pe_temp[27]={0.0}; //for the six nodes and the bubble + double Pe_gaussian[27]={0.0}; //for the six nodes and the bubble + double Ke_temp[27][3]={0.0}; //for the six nodes and the bubble + double Pe_reduced[numdof]; //for the six nodes only + double Ke_gaussian[27][3]; + double l1l6[6]; //for the six nodes of the penta + double l1l7[7]; //for the six nodes and the bubble + double B[8][27]; + double B_prime[8][27]; + double B_prime_bubble[8][3]; + double Jdet; + double Jdet2d; + double D[8][8]={0.0}; + double D_scalar; + double tBD[27][8]; + Tria* tria=NULL; + double stokesreconditioning; + int approximation; + + /*Initialize Element vector and return if necessary*/ + if(IsOnWater() || !IsOnBed() || !IsOnShelf()) return NULL; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation!=StokesApproximationEnum) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); + Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input); + + + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + bed_input->GetParameterValue(&bed, gauss); + GetNodalFunctionsP1(l1l6, gauss); + + water_pressure=gravity*rho_water*bed; + + BedNormal(&bed_normal[0],xyz_list_tria); + + for(i=0;ivalues[i*NDOF4+j]+=water_pressure*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j]; + } + + /*Clean up and return*/ + delete gauss; + return pe; } /*}}}*/ /*FUNCTION Penta::CreatePVectorAdjointStokes{{{1*/ -void Penta::CreatePVectorAdjointStokes(Vec p_g){ - - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater() || !IsOnSurface())return; - - /*Call Tria's function*/ - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->CreatePVectorAdjointStokes(p_g); +ElementVector* Penta::CreatePVectorAdjointStokes(void){ + + if (!IsOnSurface() || IsOnWater()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + ElementVector* pe=tria->CreatePVectorAdjointStokes(); delete tria->matice; delete tria; - return; + + /*clean up and return*/ + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5925) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5926) @@ -152,14 +152,20 @@ void CreatePVectorBalancedthickness( Vec pg); void CreatePVectorBalancedvelocities( Vec pg); - void CreatePVectorAdjointHoriz( Vec pg); - void CreatePVectorAdjointMacAyeal( Vec pg); - void CreatePVectorAdjointPattyn( Vec pg); - void CreatePVectorCouplingPattynStokes( Vec pg); - void CreatePVectorDiagnosticHoriz( Vec pg); + ElementVector* CreatePVectorAdjointHoriz(void); + ElementVector* CreatePVectorAdjointMacAyeal(void); + ElementVector* CreatePVectorAdjointPattyn(void); + ElementVector* CreatePVectorCouplingPattynStokes(void); + ElementVector* CreatePVectorCouplingPattynStokesViscous(void); + ElementVector* CreatePVectorCouplingPattynStokesFriction(void); + ElementVector* CreatePVectorDiagnosticHoriz(void); void CreatePVectorDiagnosticHutter( Vec pg); - void CreatePVectorDiagnosticMacAyeal( Vec pg); - void CreatePVectorDiagnosticPattyn( Vec pg); - void CreatePVectorDiagnosticStokes( Vec pg); - void CreatePVectorAdjointStokes( Vec pg); + ElementVector* CreatePVectorDiagnosticMacAyeal(void); + ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void); + ElementVector* CreatePVectorDiagnosticPattyn(void); + ElementVector* CreatePVectorDiagnosticPattynStokes(void); + ElementVector* CreatePVectorDiagnosticStokes(void); + ElementVector* CreatePVectorDiagnosticStokesViscous(void); + ElementVector* CreatePVectorDiagnosticStokesShelf(void); + ElementVector* CreatePVectorAdjointStokes(void); void CreatePVectorDiagnosticVert( Vec pg); void CreatePVectorMelting( Vec pg); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5925) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5926) @@ -3968,14 +3968,10 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/ -void Tria::CreatePVectorAdjointStokes(Vec p_g){ - +ElementVector* Tria::CreatePVectorAdjointStokes(void){ + + const int numdof=NDOF4*NUMVERTICES; int i; - - /* node data: */ - const int numdof=NDOF4*NUMVERTICES; + int ig; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; - - /* grid data: */ double vx_list[NUMVERTICES]; double vy_list[NUMVERTICES]; @@ -3985,24 +3981,11 @@ double duy_list[NUMVERTICES]; double weights_list[NUMVERTICES]; - - /* gaussian points: */ - int ig; GaussTria* gauss=NULL; - - /* parameters: */ double obs_velocity_mag,velocity_mag; double dux,duy; double meanvel, epsvel; - - /*element vector : */ double pe_g[numdof]={0.0}; - - /* Jacobian: */ double Jdet; - - /*nodal functions: */ double l1l2l3[3]; - - /*relative and algorithmic fitting: */ double scalex=0; double scaley=0; @@ -4012,13 +3995,12 @@ int response; - /*retrieve some parameters: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); this->parameters->FindParam(&meanvel,MeanVelEnum); this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /* Recover input data: */ GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); @@ -4026,5 +4008,4 @@ GetParameterListOnVertices(&vy_list[0],VyEnum); GetParameterListOnVertices(&weights_list[0],WeightsEnum); - inputs->GetParameterValue(&response,CmResponseEnum); if(response==SurfaceAverageVelMisfitEnum){ @@ -4156,29 +4137,19 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions value at gaussian point:*/ GetNodalFunctions(l1l2l3, gauss); - /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */ - - /*Compute absolute(x/y) at gaussian point: */ TriaRef::GetParameterValue(&dux, &dux_list[0],gauss); TriaRef::GetParameterValue(&duy, &duy_list[0],gauss); - /*compute Du*/ for (i=0;iweight*l1l2l3[i]; - pe_g[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; + pe->values[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i]; + pe->values[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i]; } } - - /*Add pe_g to global vector p_g: */ - VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES); /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5925) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5926) @@ -145,5 +145,5 @@ ElementVector* CreatePVectorDiagnosticMacAyeal(void); ElementVector* CreatePVectorAdjointHoriz(void); - void CreatePVectorAdjointStokes(Vec pg); + ElementVector* CreatePVectorAdjointStokes(void); ElementVector* CreatePVectorAdjointBalancedthickness(void); ElementVector* CreatePVectorDiagnosticHutter(void); Index: /issm/trunk/src/c/objects/Numerics/ElementVector.cpp =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementVector.cpp (revision 5925) +++ /issm/trunk/src/c/objects/Numerics/ElementVector.cpp (revision 5926) @@ -31,4 +31,104 @@ } /*}}}*/ +/*FUNCTION ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2){{{1*/ +ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2){ + + /*intermediaries*/ + int i,j,counter; + int gsize,fsize,ssize; + int* P=NULL; + bool found; + + /*If one of the two matrix is NULL, we copy the other one*/ + if(!pe1 && !pe2){ + ISSMERROR("Two input element matrices are NULL"); + } + else if(!pe1){ + this->Init(pe2); + return; + } + else if(!pe2){ + this->Init(pe1); + return; + } + + /*Initialize itransformation matrix pe[P[i]] = pe2[i]*/ + P=(int*)xmalloc(pe2->nrows*sizeof(int)); + + /*1: Get the new numbering of pe2 and get size of the new matrix*/ + gsize=pe1->nrows; + for(i=0;inrows;i++){ + found=false; + for(j=0;jnrows;j++){ + if(pe2->gglobaldoflist[i]==pe1->gglobaldoflist[j]){ + found=true; P[i]=j; break; + } + } + if(!found){ + P[i]=gsize; gsize++; + } + } + + /*2: Initialize static fields*/ + this->nrows=gsize; + this->pf=pe1->pf; + + /*Gset and values*/ + this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int)); + this->values=(double*)xcalloc(this->nrows,sizeof(double)); + for(i=0;inrows;i++){ + this->values[i] += pe1->values[i]; + this->gglobaldoflist[i]=pe1->gglobaldoflist[i]; + } + for(i=0;inrows;i++){ + this->values[P[i]] += pe2->values[i]; + this->gglobaldoflist[P[i]]=pe2->gglobaldoflist[i]; + } + + /*Fset*/ + fsize=pe1->fsize; + for(i=0;ifsize;i++){ + if(P[pe2->flocaldoflist[i]] >= pe1->nrows) fsize++; + } + this->fsize=fsize; + if(fsize){ + this->flocaldoflist =(int*)xmalloc(fsize*sizeof(int)); + this->fglobaldoflist=(int*)xmalloc(fsize*sizeof(int)); + for(i=0;ifsize;i++){ + this->flocaldoflist[i] =pe1->flocaldoflist[i]; + this->fglobaldoflist[i]=pe1->fglobaldoflist[i]; + } + counter=pe1->fsize; + for(i=0;ifsize;i++){ + if(P[pe2->flocaldoflist[i]] >= pe1->nrows){ + this->flocaldoflist[counter] =P[pe2->flocaldoflist[i]]; + this->fglobaldoflist[counter]=pe2->fglobaldoflist[i]; + } + } + } + else{ + this->flocaldoflist=NULL; + this->fglobaldoflist=NULL; + } + + /*clean-up*/ + xfree((void**)&P); +} +/*}}}*/ +/*FUNCTION ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2,ElementVector* pe3){{{1*/ +ElementVector::ElementVector(ElementVector* pe1, ElementVector* pe2,ElementVector* pe3){ + + /*Concatenate all matrices*/ + ElementVector* pe12 =new ElementVector(pe1,pe2); + ElementVector* pe123=new ElementVector(pe12,pe3); + + /*Initialize current object with this matrix*/ + this->Init(pe123); + + /*clean-up*/ + delete pe12; + delete pe123; +} +/*}}}*/ /*FUNCTION ElementVector::ElementVector(int gsize,int* gglobaldoflist){{{1*/ ElementVector::ElementVector(int gsize,int* in_gglobaldoflist){ @@ -49,4 +149,5 @@ } /*not needed: */ + this->fsize=0; this->flocaldoflist=NULL; this->fglobaldoflist=NULL; @@ -130,2 +231,52 @@ } /*}}}*/ +/*FUNCTION ElementVector::Echo{{{1*/ +void ElementVector::Echo(void){ + + int i,j; + printf("Element Vector echo: \n"); + printf(" nrows: %i\n",nrows); + printf(" pf: %s\n",pf?"true":"false"); + + printf(" values: \n"); + for(i=0;inrows =pe->nrows; + this->pf =pe->pf; + + this->values=(double*)xmalloc(this->nrows*sizeof(double)); + memcpy(this->values,pe->values,this->nrows*sizeof(double)); + + this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int)); + memcpy(this->gglobaldoflist,pe->gglobaldoflist,this->nrows*sizeof(int)); + + this->fsize=pe->fsize; + if(this->fsize){ + this->flocaldoflist=(int*)xmalloc(this->fsize*sizeof(int)); + memcpy(this->flocaldoflist,pe->flocaldoflist,this->fsize*sizeof(int)); + this->fglobaldoflist=(int*)xmalloc(this->fsize*sizeof(int)); + memcpy(this->fglobaldoflist,pe->fglobaldoflist,this->fsize*sizeof(int)); + } + else{ + this->flocaldoflist=NULL; + this->fglobaldoflist=NULL; + } +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Numerics/ElementVector.h =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementVector.h (revision 5925) +++ /issm/trunk/src/c/objects/Numerics/ElementVector.h (revision 5926) @@ -33,4 +33,6 @@ /*ElementVector constructors, destructors {{{1*/ ElementVector(); + ElementVector(ElementVector* pe1,ElementVector* pe2); + ElementVector(ElementVector* pe1,ElementVector* pe2,ElementVector* pe3); ElementVector(int gsize,int* gglobaldoflist); ElementVector(int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize); @@ -40,4 +42,6 @@ void AddValues(double* pe_gg); void AddToGlobal(Vec pg, Vec pf); + void Echo(void); + void Init(ElementVector* pe); /*}}}*/ };