Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 2712) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 2713) @@ -16,4 +16,5 @@ #include "../include/typedefs.h" +/*Object constructors and destructor*/ /*FUNCTION Penta constructor {{{1*/ Penta::Penta(){ @@ -70,75 +71,7 @@ } /*}}}*/ -/*FUNCTION Penta Echo {{{1*/ -void Penta::Echo(void){ - - printf("Penta:\n"); - printf(" id: %i\n",id); - printf(" mid: %i\n",mid); - printf(" mparid: %i\n",mparid); - printf(" numparid: %i\n",numparid); - - printf(" node_ids=[%i,%i,%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3],node_ids[4],node_ids[5]); - printf(" node_offsets=[%i,%i,%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3],node_offsets[4],node_offsets[5]); - printf(" matice_offset=%i\n",matice_offset); - printf(" matpar_offset=%i\n",matpar_offset); - - printf(" h=[%g,%g,%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3],h[4],h[5]); - printf(" s=[%g,%g,%g,%g,%g,%g]\n",s[0],s[1],s[2],s[3],s[4],s[5]); - printf(" b=[%g,%g,%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3],b[4],b[5]); - printf(" k=[%g,%g,%g,%g,%g,%g]\n",k[0],k[1],k[2],k[3],k[4],k[5]); - - printf(" friction_type: %i\n",friction_type); - printf(" p: %g\n",p); - printf(" q: %g\n",q); - printf(" shelf: %i\n",shelf); - printf(" onbed: %i\n",onbed); - printf(" onwater: %i\n",onwater); - printf(" onsurface: %i\n",onsurface); - printf(" collapse: %i\n",collapse); - - printf(" melting=[%g,%g,%g,%g,%g,%g]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]); - printf(" accumulation=[%g,%g,%g,%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]); - printf(" geothermalflux=[%g,%g,%g,%g,%g,%g]\n",geothermalflux[0],geothermalflux[1],geothermalflux[2],geothermalflux[3],geothermalflux[4],geothermalflux[5]); - printf(" thermal_steadystate: %i\n",thermal_steadystate); - return; -} -/*}}}*/ -/*FUNCTION Penta DeepEcho {{{1*/ -void Penta::DeepEcho(void){ - - printf("Penta:\n"); - printf(" id: %i\n",id); - printf(" mid: %i\n",mid); - printf(" mparid: %i\n",mparid); - printf(" numparid: %i\n",numparid); - - printf(" node_ids=[%i,%i,%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3],node_ids[4],node_ids[5]); - printf(" node_offsets=[%i,%i,%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3],node_offsets[4],node_offsets[5]); - printf(" matice_offset=%i\n",matice_offset); - printf(" matpar_offset=%i\n",matpar_offset); - - printf(" h=[%i,%i,%i,%i,%i,%i]\n",h[0],h[1],h[2],h[3],h[4],h[5]); - printf(" s=[%i,%i,%i,%i,%i,%i]\n",s[0],s[1],s[2],s[3],s[4],s[5]); - printf(" b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]); - printf(" k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]); - - printf(" friction_type: %i\n",friction_type); - printf(" p: %g\n",p); - printf(" q: %g\n",q); - printf(" shelf: %i\n",shelf); - printf(" onbed: %i\n",onbed); - printf(" onwater: %i\n",onwater); - printf(" onsurface: %i\n",onsurface); - printf(" collapse: %i\n",collapse); - - printf(" melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]); - printf(" accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]); - printf(" geothermalflux=[%i,%i,%i,%i,%i,%i]\n",geothermalflux[0],geothermalflux[1],geothermalflux[2],geothermalflux[3],geothermalflux[4],geothermalflux[5]); - printf(" thermal_steadystate: %i\n",thermal_steadystate); - return; -} -/*}}}*/ -/*FUNCTION Penta Marshall {{{1*/ + +/*Object marshall*/ +/*FUNCTION Marshall {{{1*/ void Penta::Marshall(char** pmarshalled_dataset){ @@ -190,5 +123,5 @@ } /*}}}*/ -/*FUNCTION Penta MarshallSize {{{1*/ +/*FUNCTION MarshallSize {{{1*/ int Penta::MarshallSize(){ @@ -225,10 +158,5 @@ } /*}}}*/ -/*FUNCTION Penta GetName {{{1*/ -char* Penta::GetName(void){ - return "penta"; -} -/*}}}*/ -/*FUNCTION Penta Demarshall {{{1*/ +/*FUNCTION Demarshall {{{1*/ void Penta::Demarshall(char** pmarshalled_dataset){ @@ -283,24 +211,44 @@ } /*}}}*/ -/*FUNCTION Penta Enum {{{1*/ -int Penta::Enum(void){ - - return PentaEnum(); - -} -/*}}}*/ -/*FUNCTION Penta GetId {{{1*/ -int Penta::GetId(void){ - return id; -} -/*}}}*/ -/*FUNCTION Penta MyRank {{{1*/ -int Penta::MyRank(void){ - extern int my_rank; - return my_rank; -} -/*}}}*/ -/*FUNCTION Penta Configure {{{1*/ -#undef __FUConfigure NCT__ + +/*Object functions*/ +/*FUNCTION ComputePressure {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::ComputePressure" +void Penta::ComputePressure(Vec pg){ + + int i; + const int numgrids=6; + int doflist[numgrids]; + double pressure[numgrids]; + double rho_ice,g; + double xyz_list[numgrids][3]; + + /*If on water, skip: */ + if(onwater)return; + + /*Get node data: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + + /*pressure is lithostatic: */ + //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab + + /*Get dof list on which we will plug the pressure values: */ + GetDofList1(&doflist[0]); + + /*pressure is lithostatic: */ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + for(i=0;iCreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; - } - - /*Now, onto the formulation for the vertical velocity: */ - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;iCreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + } + + /*Now, onto the formulation for the vertical velocity: */ + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iCreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION CreateKMatrixPrognostic {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixPrognostic" + +void Penta::CreateKMatrixPrognostic(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; + +} +/*}}}*/ +/*FUNCTION CreateKMatrixSlopeCompute {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixSlopeCompute" + +void Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; + +} +/*}}}*/ +/*FUNCTION CreateKMatrixThermal {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixThermal" +void Penta::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ + + /* local declarations */ + int i,j; + int found=0; + + /* node data: */ + const int numgrids=6; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_area_gauss,igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + double gauss_l1l2l3[3]; + + int area_order=5; + int num_vert_gauss=5; + + int dofs[3]={0,1,2}; + double dt; + double K[2][2]={0.0}; + + double vxvyvz_list[numgrids][3]; + double vx_list[numgrids]; + int vx_list_exists; + double u; + double vy_list[numgrids]; + int vy_list_exists; + double v; + double vz_list[numgrids]; + int vz_list_exists; + double w; + + /*matrices: */ + double K_terms[numdof][numdof]={0.0}; + double Ke_gaussian_conduct[numdof][numdof]; + double Ke_gaussian_advec[numdof][numdof]; + double Ke_gaussian_artdiff[numdof][numdof]; + double Ke_gaussian_transient[numdof][numdof]; + double B[3][numdof]; + double Bprime[3][numdof]; + double B_conduct[3][numdof]; + double B_advec[3][numdof]; + double B_artdiff[2][numdof]; + double Bprime_advec[3][numdof]; + double L[numdof]; + double D_scalar; + double D[3][3]; + double l1l2l3[3]; + double tl1l2l3D[3]; + double tBD[3][numdof]; + double tBD_conduct[3][numdof]; + double tBD_advec[3][numdof]; + double tBD_artdiff[3][numdof]; + double tLD[numdof]; + + double Jdet; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + double heatcapacity,thermalconductivity; + double mixed_layer_capacity,thermal_exchange_velocity; + + /*Collapsed formulation: */ + Tria* tria=NULL; + ParameterInputs* inputs=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + // /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + thermalconductivity=matpar->GetThermalConductivity(); + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + for(i=0;iartdiff){ + /*Build K: */ + D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+numpar->epsvel); + if(dt){ + D_scalar=D_scalar*dt; + } + K[0][0]=D_scalar*pow(u,2); K[0][1]=D_scalar*fabs(u)*fabs(v); + K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2); + + /*Get B_artdiff: */ + GetB_artdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss_coord); + + /* Do the triple product B'*K*B: */ + MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0); + MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0); + } + else{ + for(i=0;iCreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + } +} +/*}}}*/ +/*FUNCTION CreatePVector {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVector" @@ -1116,719 +1393,5 @@ } /*}}}*/ -/*FUNCTION Penta UpdateFromInputs {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::UpdateFromInputs" -void Penta::UpdateFromInputs(void* vinputs){ - - int dofs[1]={0}; - double temperature_list[6]; - double temperature_average; - double B_list[6]; - double B_average; - - ParameterInputs* inputs=NULL; - - /*If on water, skip: */ - if(onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*Update internal data if inputs holds new values: */ - inputs->Recover("thickness",&h[0],1,dofs,6,(void**)nodes); - inputs->Recover("surface",&s[0],1,dofs,6,(void**)nodes); - inputs->Recover("bed",&b[0],1,dofs,6,(void**)nodes); - inputs->Recover("drag",&k[0],1,dofs,6,(void**)nodes); - inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes); - inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes); - - //Update material if necessary - if(inputs->Recover("temperature",&temperature_list[0],1,dofs,6,(void**)nodes)){ - if(matice && !collapse){ - //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2]) - // +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0; - temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; - B_average=Paterson(temperature_average); - matice->SetB(B_average); - } - } - - if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){ - if(matice && collapse){ - temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; - B_average=Paterson(temperature_average); - //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2]) - // +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0; - matice->SetB(B_average); - } - } - - if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){ - if(matice){ - B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0; - matice->SetB(B_average); - } - } - -} -/*}}}*/ -/*FUNCTION Penta GetMatPar {{{1*/ -void* Penta::GetMatPar(){ - return matpar; -} -/*}}}*/ -/*FUNCTION Penta GetShelf {{{1*/ -int Penta::GetShelf(){ - return shelf; -} -/*}}}*/ -/*FUNCTION Penta GetNodes {{{1*/ -void Penta::GetNodes(void** vpnodes){ - int i; - Node** pnodes=(Node**)vpnodes; - - for(i=0;i<6;i++){ - pnodes[i]=nodes[i]; - } -} -/*}}}*/ -/*FUNCTION Penta GetOnBed {{{1*/ -int Penta::GetOnBed(){ - return onbed; -} -/*}}}*/ -/*FUNCTION Penta GetThicknessList {{{1*/ -void Penta::GetThicknessList(double* thickness_list){ - - int i; - for(i=0;i<6;i++)thickness_list[i]=h[i]; -} -/*}}}*/ -/*FUNCTION Penta GetBedList {{{1*/ -void Penta::GetBedList(double* bed_list){ - - int i; - for(i=0;i<6;i++)bed_list[i]=b[i]; - -} -/*}}}*/ -/*FUNCTION Penta copy {{{1*/ -Object* Penta::copy() { - return new Penta(*this); -} -/*}}}*/ -/*FUNCTION Penta Du {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::Du" -void Penta::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){ - - int i; - Tria* tria=NULL; - - /*If on water, skip: */ - if(onwater)return; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return; - } - else if (collapse){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute Du*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->Du(du_g,inputs,analysis_type,sub_analysis_type); - delete tria; - return; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->Du(du_g,inputs,analysis_type,sub_analysis_type); - delete tria; - return; - } -} -/*}}}*/ -/*FUNCTION Penta Gradj {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::Gradj" -void Penta::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ - - /*If on water, skip grad (=0): */ - if(onwater)return; - - if (strcmp(control_type,"drag")==0){ - GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type); - } - else if (strcmp(control_type,"B")==0){ - GradjB( grad_g, inputs,analysis_type,sub_analysis_type); - } - else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type)); -} -/*}}}*/ -/*FUNCTION Penta GradjDrag {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::GradjDrag" -void Penta::GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ - - Tria* tria=NULL; - - /*If on water, skip: */ - if(onwater)return; - - /*If on shelf, skip: */ - if(shelf)return; - - /*Bail out if this element does not touch the bedrock: */ - if (!onbed) return; - - if (sub_analysis_type==HorizAnalysisEnum()){ - - /*MacAyeal or Pattyn*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type); - delete tria; - return; - } - else if (sub_analysis_type==StokesAnalysisEnum()){ - - /*Stokes*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type); - delete tria; - return; - } - else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); -} -/*}}}*/ -/*FUNCTION Penta GradjB {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::GradjB" -void Penta::GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ - - Tria* tria=NULL; - - /*If on water, skip: */ - if(onwater)return; - - if (collapse){ - /*Bail out element if collapsed (2d) and not on bed*/ - if (!onbed) return; - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute gardj*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type); - delete tria; - return; - } - else{ - /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type); - delete tria; - return; - } -} -/*}}}*/ -/*FUNCTION Penta Misfit {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::Misfit" -double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){ - - double J; - Tria* tria=NULL; - - /*If on water, return 0: */ - if(onwater)return 0; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return 0; - } - else if (collapse){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute Misfit*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->Misfit(inputs,analysis_type,sub_analysis_type); - delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->Misfit(inputs,analysis_type,sub_analysis_type); - delete tria; - return J; - } -} -/*}}}*/ -/*FUNCTION Penta SpawnTria {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::SpawnTria" -void* Penta::SpawnTria(int g0, int g1, int g2){ - - /*out of grids g0,g1 and g2 from Penta, build a tria element: */ - Tria* tria=NULL; - double tria_h[3]; - double tria_s[3]; - double tria_b[3]; - double tria_c[3]; - double tria_k[3]; - double tria_melting[3]; - double tria_accumulation[3]; - double tria_geothermalflux[3]; - - /*configuration: */ - int tria_node_ids[3]; - Node* tria_nodes[3]; - int tria_node_offsets[3]; - - tria_h[0]=h[g0]; - tria_h[1]=h[g1]; - tria_h[2]=h[g2]; - - tria_s[0]=s[g0]; - tria_s[1]=s[g1]; - tria_s[2]=s[g2]; - - tria_b[0]=b[g0]; - tria_b[1]=b[g1]; - tria_b[2]=b[g2]; - - tria_k[0]=k[g0]; - tria_k[1]=k[g1]; - tria_k[2]=k[g2]; - - tria_melting[0]=melting[g0]; - tria_melting[1]=melting[g1]; - tria_melting[2]=melting[g2]; - - tria_accumulation[0]=accumulation[g0]; - tria_accumulation[1]=accumulation[g1]; - tria_accumulation[2]=accumulation[g2]; - - tria_geothermalflux[0]=geothermalflux[g0]; - tria_geothermalflux[1]=geothermalflux[g1]; - tria_geothermalflux[2]=geothermalflux[g2]; - - tria_node_ids[0]=node_ids[g0]; - tria_node_ids[1]=node_ids[g1]; - tria_node_ids[2]=node_ids[g2]; - - tria_nodes[0]=nodes[g0]; - tria_nodes[1]=nodes[g1]; - tria_nodes[2]=nodes[g2]; - - tria_node_offsets[0]=node_offsets[g0]; - tria_node_offsets[1]=node_offsets[g1]; - tria_node_offsets[2]=node_offsets[g2]; - - tria= new Tria(id,mid,mparid,numparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,onwater); - - tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets); - tria->MaticeConfiguration(matice,matice_offset); - tria->MatparConfiguration(matpar,matpar_offset); - tria->NumparConfiguration(numpar,numpar_offset); - - return tria; - -} -/*}}}*/ -/*FUNCTION Penta GetDofList {{{1*/ -void Penta::GetDofList(int* doflist,int* pnumberofdofspernode){ - - int i,j; - int doflist_per_node[MAXDOFSPERNODE]; - int numberofdofspernode; - - for(i=0;i<6;i++){ - nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); - for(j=0;jGetDofList1(); - } - -} -/*}}}*/ -/*FUNCTION Penta GetStrainRate {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::GetStrainRate" -void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4){ - - int i; - const int numgrids=6; - const int NDOF2=2; - - double B[5][NDOF2*numgrids]; - - /*Get B matrix: */ - GetB(&B[0][0], xyz_list, gauss_l1l2l3l4); - -#ifdef _ISSM_DEBUG_ - printf("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); - printf(" [ %lf %lf \n",B[1][0],B[1][1]); - printf(" [ %lf %lf ]\n",B[2][0],B[2][1]); - printf(" [ %lf %lf ]\n",B[3][0],B[3][1]); - printf(" [ %lf %lf ]\n",B[4][0],B[4][1]); - - printf("B for grid2 : [ %lf %lf \n",B[0][2],B[0][3]); - printf(" [ %lf %lf \n",B[1][2],B[1][3]); - printf(" [ %lf %lf ]\n",B[2][2],B[2][3]); - printf(" [ %lf %lf ]\n",B[3][2],B[3][3]); - printf(" [ %lf %lf ]\n",B[4][2],B[4][3]); - - printf("B for grid3 : [ %lf %lf \n", B[0][4],B[0][5]); - printf(" [ %lf %lf \n", B[1][4],B[1][5]); - printf(" [ %lf %lf ]\n",B[2][4],B[2][5]); - printf(" [ %lf %lf ]\n",B[3][4],B[3][5]); - printf(" [ %lf %lf ]\n",B[4][4],B[4][5]); - - printf("B for grid4 : [ %lf %lf \n", B[0][6],B[0][7]); - printf(" [ %lf %lf \n", B[1][6],B[1][7]); - printf(" [ %lf %lf ]\n",B[2][6],B[2][7]); - printf(" [ %lf %lf ]\n",B[3][6],B[3][7]); - printf(" [ %lf %lf ]\n",B[4][6],B[4][7]); - - printf("B for grid5 : [ %lf %lf \n", B[0][8],B[0][9]); - printf(" [ %lf %lf \n", B[1][8],B[1][9]); - printf(" [ %lf %lf ]\n",B[2][8],B[2][9]); - printf(" [ %lf %lf ]\n",B[3][8],B[3][9]); - printf(" [ %lf %lf ]\n",B[4][8],B[4][9]); - - printf("B for grid6 : [ %lf %lf \n", B[0][10],B[0][11]); - printf(" [ %lf %lf \n", B[1][10],B[1][11]); - printf(" [ %lf %lf ]\n",B[2][10],B[2][11]); - printf(" [ %lf %lf ]\n",B[3][10],B[3][11]); - printf(" [ %lf %lf ]\n",B[4][10],B[4][11]); - -#endif - - /*Multiply B by velocity, to get strain rate: */ - MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, - velocity,NDOF2*numgrids,1,0, - epsilon,0); - -} -/*}}}*/ -/*FUNCTION Penta GetB {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::GetB" -void Penta::GetB(double* B, double* xyz_list, double* gauss_l1l2l3l4){ - - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. - * For grid i, Bi can be expressed in the basic coordinate system - * by: - * Bi_basic=[ dh/dx 0 ] - * [ 0 dh/dy ] - * [ 1/2*dh/dy 1/2*dh/dx ] - * [ 1/2*dh/dz 0 ] - * [ 0 1/2*dh/dz ] - * where h is the interpolation function for grid i. - * - * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) - */ - - int i; - const int numgrids=6; - const int NDOF3=3; - const int NDOF2=2; - - double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids]; - - /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */ - GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); - -#ifdef _ISSM_DEBUG_ - for (i=0;iGetDofList(&nodedofs[0],&numberofdofspernode); - VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES); - VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); + const int DOFPERGRID=4; + const int numdof=numgrids*DOFPERGRID; + const int numgrids2d=3; + int numdof2d=numgrids2d*DOFPERGRID; + int doflist[numdof]; + int numberofdofspernode; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + + /*parameters: */ + double xyz_list_tria[numgrids2d][3]; + double xyz_list[numgrids][3]; + double surface_normal[3]; + double bed_normal[3]; + double bed; + double vxvyvz_list[numgrids][3]; + + /* gaussian points: */ + int num_area_gauss; + int igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + + /* specific gaussian point: */ + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + double gauss_coord_tria[3]; + + int area_order=5; + int num_vert_gauss=5; + + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double viscosity; + double water_pressure; + int dofs[3]={0,1,2}; + + /*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 + 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 L[3]; //for the three nodes of the bed + 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]; + double P_terms[numdof]; + + ParameterInputs* inputs=NULL; + Tria* tria=NULL; + + /*If on water, skip load: */ + if(onwater)return; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Set P_terms to 0: */ + for(i=0;iGetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + + /*recover extra inputs from users, at current convergence iteration: */ + inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetViscosity3dStokes(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /* Get nodal functions */ + GetNodalFunctionsStokes(&l1l7[0], gauss_coord); + + /* Build gaussian vector */ + for(i=0;iGetDofList1(); - - /*get field for this base node: */ - fieldel=field_serial[dof1]; - - //go throfieldn all nodes which sit on top of this node, until we reach the surface, - //and plfield field in field - for(;;){ - - dof1=node->GetDofList1(); - VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + * onto this scalar matrix, so that we win some computational time: */ + D_scalar=gauss_weight*Jdet; + for (i=0;i<6;i++){ + D[i][i]=D_scalar*viscosity; + } + for (i=6;i<8;i++){ + D[i][i]=-D_scalar*numpar->stokesreconditioning; + } + + /* Do the triple product tB*D*Bprime: */ + MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); + MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0); + + /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ + for(i=0;i<27;i++){ + for(j=0;j<3;j++){ + Ke_temp[i][j]+=Ke_gaussian[i][j]; } } } - else if ( - (strcmp(field_name,"thickness")==0) || - (strcmp(field_name,"surface")==0) || - (strcmp(field_name,"bed")==0) || - (strcmp(field_name,"slopex")==0) || - (strcmp(field_name,"slopey")==0) - ){ - - /* node data: */ - const int numdof=1*numgrids; - int doflist[numdof]; - int nodedofs; - double fieldel; - - GetDofList(&doflist[0],&numberofdofspernode); - - /*this penta is on the bed. For each node on the base of this penta, - * we grab the thickness. Once we know the thickness, we follow the upper nodes, - * inserting the same thickness value into tg, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get velocity for this base node: */ - fieldel=field_serial[doflist[numberofdofspernode*i+0]]; - - //go through all nodes which sit on top of this node, until we reach the surface, - //and pltg fieldel in field: - for(;;){ - - node->GetDofList(&nodedofs,&numberofdofspernode); - VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); + } + + /*Deal with 2d friction at the bedrock interface: */ + if ( (onbed==1) && (shelf==1)){ + + for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); + + /* Get bed at gaussian point */ + GetParameterValue(&bed,&b[0],gauss_coord); + + /*Get L matrix : */ + tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1); + + /*Get water_pressure at gaussian point */ + water_pressure=gravity*rho_water*bed; + + /*Get normal vecyor to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1]=-surface_normal[1]; + bed_normal[2]=-surface_normal[2]; + + for(i=0;iGetRhoIce(); - g=matpar->GetG(); - for(i=0;iCreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; +/*FUNCTION CreatePVectorMelting {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorMelting" +void Penta::CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ return; - -} -/*}}}*/ -/*FUNCTION Penta CreatePVectorSlopeCompute {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::CreatePVectorSlopeCompute" - -void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +} +/*}}}*/ +/*FUNCTION CreatePVectorPrognostic {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorPrognostic" + +void Penta::CreatePVectorPrognostic( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -2462,32 +1967,9 @@ } /*}}}*/ -/*FUNCTION Penta CreateKMatrixPrognostic {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::CreateKMatrixPrognostic" - -void Penta::CreateKMatrixPrognostic(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ - - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(onwater)return; - - /*Is this element on the bed? :*/ - if(!onbed)return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; - return; - -} -/*}}}*/ -/*FUNCTION Penta CreatePVectorPrognostic {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::CreatePVectorPrognostic" - -void Penta::CreatePVectorPrognostic( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +/*FUNCTION CreatePVectorSlopeCompute {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorSlopeCompute" + +void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -2507,156 +1989,833 @@ } /*}}}*/ -/*FUNCTION Penta ReduceMatrixStokes {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "ReduceMatrixStokes" -void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){ - +/*FUNCTION CreatePVectorThermal {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorThermal" +void Penta::CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ + + + /*indexing: */ int i,j; - - double Kii[24][24]; - double Kib[24][3]; - double Kbb[3][3]; - double Kbi[3][24]; - double Kbbinv[3][3]; - double KibKbbinv[24][3]; - double Kright[24][24]; - - /*Create the four matrices used for reduction */ - for(i=0;i<24;i++){ - for(j=0;j<24;j++){ - Kii[i][j]=*(Ke_temp+27*i+j); - } - for(j=0;j<3;j++){ - Kib[i][j]=*(Ke_temp+27*i+j+24); - } - } - for(i=0;i<3;i++){ - for(j=0;j<24;j++){ - Kbi[i][j]=*(Ke_temp+27*(i+24)+j); - } - for(j=0;j<3;j++){ - Kbb[i][j]=*(Ke_temp+27*(i+24)+j+24); - } - } - - /*Inverse the matrix corresponding to bubble part Kbb */ - GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]); - - /*Multiply matrices to create the reduce matrix Ke_reduced */ - MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0); - MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Kbi[0][0],3,24,0,&Kright[0][0],0); - - /*Affect value to the reduced matrix */ - for(i=0;i<24;i++){ - for(j=0;j<24;j++){ - *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j]; - } - } -} -/*}}}*/ -/*FUNCTION Penta GetMatrixInvert {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "GetMatrixInvert" -void Penta::GetMatrixInvert(double* Ke_invert, double* Ke){ - /*Inverse a 3 by 3 matrix only */ - - double a,b,c,d,e,f,g,h,i; - double det; + int found=0; + + const int numgrids=6; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + int doflist[numdof]; + int numberofdofspernode; + + /*Grid data: */ + double xyz_list[numgrids][3]; + + /* gaussian points: */ + int num_area_gauss,igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + int area_order=2; + int num_vert_gauss=3; + + double dt; + double vx_list[numgrids]; + double vy_list[numgrids]; + double vz_list[numgrids]; + double vxvyvz_list[numgrids][3]; + double temperature_list[numgrids]; + double temperature; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + double mixed_layer_capacity,heatcapacity; + double beta,meltingpoint,thermal_exchange_velocity; + + /* element parameters: */ + int friction_type; + + int dofs[3]={0,1,2}; + int dofs1[1]={0}; + + /*matrices: */ + double P_terms[numdof]={0.0}; + double L[numdof]; + double l1l2l3[3]; + double alpha2_list[3]; + double basalfriction_list[3]={0.0}; + double basalfriction; + double epsilon[6]; + double epsilon_sqr[3][3]; + double epsilon_matrix[3][3]; + + double Jdet; + double viscosity; + double epsilon_eff; + double phi; + double t_pmp; + double scalar; + double scalar_def; + double scalar_ocean; + double scalar_transient; + + /*Collapsed formulation: */ + Tria* tria=NULL; + ParameterInputs* inputs=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + beta=matpar->GetBeta(); + meltingpoint=matpar->GetMeltingPoint(); + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + if(dt){ + found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); + } + + for(i=0;iGetViscosity3dStokes(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /* Get nodal functions */ + GetNodalFunctions(&L[0], gauss_coord); + + /*Build deformational heating: */ + GetPhi(&phi, &epsilon[0], viscosity); + + /*Build pe_gaussian */ + scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; + if(dt){ + scalar_def=scalar_def*dt; + } + + for(i=0;iCreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type); + delete tria; + } + + /* Geothermal flux on ice sheet base and basal friction */ + if(onbed && !shelf){ + + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type); + delete tria; + } + extern int my_rank; + +cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&vert_gauss_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + +} +/*}}}*/ +/*FUNCTION DeepEcho {{{1*/ +void Penta::DeepEcho(void){ + + printf("Penta:\n"); + printf(" id: %i\n",id); + printf(" mid: %i\n",mid); + printf(" mparid: %i\n",mparid); + printf(" numparid: %i\n",numparid); + + printf(" node_ids=[%i,%i,%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3],node_ids[4],node_ids[5]); + printf(" node_offsets=[%i,%i,%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3],node_offsets[4],node_offsets[5]); + printf(" matice_offset=%i\n",matice_offset); + printf(" matpar_offset=%i\n",matpar_offset); + + printf(" h=[%i,%i,%i,%i,%i,%i]\n",h[0],h[1],h[2],h[3],h[4],h[5]); + printf(" s=[%i,%i,%i,%i,%i,%i]\n",s[0],s[1],s[2],s[3],s[4],s[5]); + printf(" b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]); + printf(" k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]); + + printf(" friction_type: %i\n",friction_type); + printf(" p: %g\n",p); + printf(" q: %g\n",q); + printf(" shelf: %i\n",shelf); + printf(" onbed: %i\n",onbed); + printf(" onwater: %i\n",onwater); + printf(" onsurface: %i\n",onsurface); + printf(" collapse: %i\n",collapse); + + printf(" melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]); + printf(" accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]); + printf(" geothermalflux=[%i,%i,%i,%i,%i,%i]\n",geothermalflux[0],geothermalflux[1],geothermalflux[2],geothermalflux[3],geothermalflux[4],geothermalflux[5]); + printf(" thermal_steadystate: %i\n",thermal_steadystate); + return; +} +/*}}}*/ +/*FUNCTION Du {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::Du" +void Penta::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){ + + int i; + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute Du*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->Du(du_g,inputs,analysis_type,sub_analysis_type); + delete tria; + return; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + tria->Du(du_g,inputs,analysis_type,sub_analysis_type); + delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION Echo {{{1*/ +void Penta::Echo(void){ + + printf("Penta:\n"); + printf(" id: %i\n",id); + printf(" mid: %i\n",mid); + printf(" mparid: %i\n",mparid); + printf(" numparid: %i\n",numparid); + + printf(" node_ids=[%i,%i,%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3],node_ids[4],node_ids[5]); + printf(" node_offsets=[%i,%i,%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3],node_offsets[4],node_offsets[5]); + printf(" matice_offset=%i\n",matice_offset); + printf(" matpar_offset=%i\n",matpar_offset); + + printf(" h=[%g,%g,%g,%g,%g,%g]\n",h[0],h[1],h[2],h[3],h[4],h[5]); + printf(" s=[%g,%g,%g,%g,%g,%g]\n",s[0],s[1],s[2],s[3],s[4],s[5]); + printf(" b=[%g,%g,%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3],b[4],b[5]); + printf(" k=[%g,%g,%g,%g,%g,%g]\n",k[0],k[1],k[2],k[3],k[4],k[5]); + + printf(" friction_type: %i\n",friction_type); + printf(" p: %g\n",p); + printf(" q: %g\n",q); + printf(" shelf: %i\n",shelf); + printf(" onbed: %i\n",onbed); + printf(" onwater: %i\n",onwater); + printf(" onsurface: %i\n",onsurface); + printf(" collapse: %i\n",collapse); + + printf(" melting=[%g,%g,%g,%g,%g,%g]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]); + printf(" accumulation=[%g,%g,%g,%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]); + printf(" geothermalflux=[%g,%g,%g,%g,%g,%g]\n",geothermalflux[0],geothermalflux[1],geothermalflux[2],geothermalflux[3],geothermalflux[4],geothermalflux[5]); + printf(" thermal_steadystate: %i\n",thermal_steadystate); + return; +} +/*}}}*/ +/*FUNCTION Enum {{{1*/ +int Penta::Enum(void){ + + return PentaEnum(); + +} +/*}}}*/ +/*FUNCTION FieldExtrude {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::FieldExtrude" +void Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){ + + /* node data: */ + const int numgrids=6; + int numberofdofspernode; + Node* node=NULL; + int i; + int extrude=0; + + /*Figure out if we should extrude for this element: */ + if (iscollapsed){ + /*From higher level, we are told to extrude only elements that have the collapse flag on: */ + if (collapse)extrude=1; + else extrude=0; + } + else{ + /*From higher level, we are told to extrude all elements: */ + extrude=1; + } + + /*Now, extrusion starts from the bed on, so double check this element is on + * the bedrock: */ + if(onbed==0)extrude=0; + + /*Go on and extrude field: */ + if (extrude){ + + if (strcmp(field_name,"velocity")==0){ + + /* node data: */ + const int numdof=2*numgrids; + int doflist[numdof]; + int nodedofs[2]; + double fieldel[2]; + + + GetDofList(&doflist[0],&numberofdofspernode); + + /*this penta is a collapsed macayeal. For each node on the base of this penta, + * we grab the field. Once we know the field, we follow the upper nodes, + * inserting the same field value into field, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + + /*get field for this base node: */ + fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]]; + fieldel[1]=field_serial[doflist[numberofdofspernode*i+1]]; + + //go throfieldn all nodes which sit on top of this node, until we reach the surface, + //and plfield field in field + for(;;){ + + node->GetDofList(&nodedofs[0],&numberofdofspernode); + VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES); + VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + } //if (strcmp(field_name,"velocity")==0) + else if (strcmp(field_name,"gradj")==0){ + + /* node data: */ + int dof1; + double fieldel; + + /*this penta is a collapsed macayeal. For each node on the base of this penta, + * we grab the field. Once we know the field, we follow the upper nodes, + * inserting the same field value into field, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + dof1=node->GetDofList1(); + + /*get field for this base node: */ + fieldel=field_serial[dof1]; + + //go throfieldn all nodes which sit on top of this node, until we reach the surface, + //and plfield field in field + for(;;){ + + dof1=node->GetDofList1(); + VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + } + else if ( + (strcmp(field_name,"thickness")==0) || + (strcmp(field_name,"surface")==0) || + (strcmp(field_name,"bed")==0) || + (strcmp(field_name,"slopex")==0) || + (strcmp(field_name,"slopey")==0) + ){ + + /* node data: */ + const int numdof=1*numgrids; + int doflist[numdof]; + int nodedofs; + double fieldel; + + GetDofList(&doflist[0],&numberofdofspernode); + + /*this penta is on the bed. For each node on the base of this penta, + * we grab the thickness. Once we know the thickness, we follow the upper nodes, + * inserting the same thickness value into tg, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + + /*get velocity for this base node: */ + fieldel=field_serial[doflist[numberofdofspernode*i+0]]; + + //go through all nodes which sit on top of this node, until we reach the surface, + //and pltg fieldel in field: + for(;;){ + + node->GetDofList(&nodedofs,&numberofdofspernode); + VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + + } + else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!")); + + } //if (extrude) +} +/*}}}*/ +/*FUNCTION GetB {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetB" +void Penta::GetB(double* B, double* xyz_list, double* gauss_l1l2l3l4){ + + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. + * For grid i, Bi can be expressed in the basic coordinate system + * by: + * Bi_basic=[ dh/dx 0 ] + * [ 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * [ 1/2*dh/dz 0 ] + * [ 0 1/2*dh/dz ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) + */ + + int i; + const int numgrids=6; + const int NDOF3=3; + const int NDOF2=2; + + double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids]; + + /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */ + GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); + +#ifdef _ISSM_DEBUG_ + for (i=0;iGetDofList(&doflist_per_node[0],&numberofdofspernode); + for(j=0;jGetDofList1(); + } + +} +/*}}}*/ +/*FUNCTION GetId {{{1*/ +int Penta::GetId(void){ + return id; +} +/*}}}*/ +/*FUNCTION GetJacobian {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetJacobian" +void Penta::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3l4){ + + const int NDOF3=3; + int i,j; + + /*The Jacobian is constant over the element, discard the gaussian points. + * J is assumed to have been allocated of size NDOF2xNDOF2.*/ + + double A1,A2,A3; //area coordinates + double xi,eta,zi; //parametric coordinates + + double x1,x2,x3,x4,x5,x6; + double y1,y2,y3,y4,y5,y6; + double z1,z2,z3,z4,z5,z6; + + double sqrt3=sqrt(3.0); + + /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */ + A1=gauss_l1l2l3l4[0]; + A2=gauss_l1l2l3l4[1]; + A3=gauss_l1l2l3l4[2]; + + xi=A2-A1; + eta=sqrt3*A3; + zi=gauss_l1l2l3l4[3]; + + x1=*(xyz_list+3*0+0); + x2=*(xyz_list+3*1+0); + x3=*(xyz_list+3*2+0); + x4=*(xyz_list+3*3+0); + x5=*(xyz_list+3*4+0); + x6=*(xyz_list+3*5+0); + + y1=*(xyz_list+3*0+1); + y2=*(xyz_list+3*1+1); + y3=*(xyz_list+3*2+1); + y4=*(xyz_list+3*3+1); + y5=*(xyz_list+3*4+1); + y6=*(xyz_list+3*5+1); + + z1=*(xyz_list+3*0+2); + z2=*(xyz_list+3*1+2); + z3=*(xyz_list+3*2+2); + z4=*(xyz_list+3*3+2); + z5=*(xyz_list+3*4+2); + z6=*(xyz_list+3*5+2); + + + *(J+NDOF3*0+0)=1.0/4.0*(x1-x2-x4+x5)*zi+1.0/4.0*(-x1+x2-x4+x5); + *(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); + *(J+NDOF3*2+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +1.0/4.0*(-x1+x5-x2+x4); + + *(J+NDOF3*0+1)=1.0/4.0*(y1-y2-y4+y5)*zi+1.0/4.0*(-y1+y2-y4+y5); + *(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); + *(J+NDOF3*2+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+1.0/4.0*(y1-y2-y4+y5)*xi+1.0/4.0*(y4-y1+y5-y2); + + *(J+NDOF3*0+2)=1.0/4.0*(z1-z2-z4+z5)*zi+1.0/4.0*(-z1+z2-z4+z5); + *(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); + *(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4); + +#ifdef _ISSM_DEBUG_ + for(i=0;i<3;i++){ + for (j=0;j<3;j++){ + printf("%lf ",*(J+NDOF3*i+j)); + } + printf("\n"); + } #endif - - /*B_primeuild B_prime: */ - for (i=0;iGetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore - get tria gaussian points as well as segment gaussian points. For tria gaussian - points, order of integration is 2, because we need to integrate the product tB*D*B' - which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian - points, same deal, which yields 3 gaussian points.*/ - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); - - /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /* Get nodal functions */ - GetNodalFunctionsStokes(&l1l7[0], gauss_coord); - - /* Build gaussian vector */ - for(i=0;istokesreconditioning; - } - - /* Do the triple product tB*D*Bprime: */ - MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); - MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0); - - /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ - for(i=0;i<27;i++){ - for(j=0;j<3;j++){ - Ke_temp[i][j]+=Ke_gaussian[i][j]; - } - } - } - } - - /*Deal with 2d friction at the bedrock interface: */ - if ( (onbed==1) && (shelf==1)){ - - for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); - - /* Get bed at gaussian point */ - GetParameterValue(&bed,&b[0],gauss_coord); - - /*Get L matrix : */ - tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1); - - /*Get water_pressure at gaussian point */ - water_pressure=gravity*rho_water*bed; - - /*Get normal vecyor to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list_tria); - - bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result - bed_normal[1]=-surface_normal[1]; - bed_normal[2]=-surface_normal[2]; - - for(i=0;iGetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - heatcapacity=matpar->GetHeatCapacity(); - thermalconductivity=matpar->GetThermalConductivity(); - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - for(i=0;iartdiff){ - /*Build K: */ - D_scalar=gauss_weight*Jdet/(pow(u,2)+pow(v,2)+numpar->epsvel); - if(dt){ - D_scalar=D_scalar*dt; - } - K[0][0]=D_scalar*pow(u,2); K[0][1]=D_scalar*fabs(u)*fabs(v); - K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2); - - /*Get B_artdiff: */ - GetB_artdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss_coord); - - /* Do the triple product B'*K*B: */ - MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0); - MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0); - } - else{ - for(i=0;iCreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; - } -} -/*}}}*/ -/*FUNCTION Penta GetB_conduct {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::GetB_conduct" -void Penta::GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord){ - - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. - * For grid i, Bi' can be expressed in the basic coordinate system - * by: - * Bi_conduct_basic=[ dh/dx ] - * [ dh/dy ] - * [ dh/dz ] - * where h is the interpolation function for grid i. +/*FUNCTION GetParameterDerivativeValue {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetParameterDerivativeValue" +void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){ + + /*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_l1l2l3l4: + * 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; + * 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; + * 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; * - * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) + * p is a vector of size 3x1 already allocated. */ + const int NDOF3=3; + const int numgrids=6; + double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids]; + + /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */ + GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); + + *(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5]; + ; + *(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5]; + + *(p+2)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[2][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[2][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[2][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[2][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[2][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[2][5]; + +} +/*}}}*/ +/*FUNCTION GetParameterValue {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetParameterValue" +void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){ + + const int numgrids=6; + double l1l2l3l4l5l6[numgrids]; + + GetNodalFunctions(&l1l2l3l4l5l6[0], gauss_l1l2l3l4); + + *pvalue=l1l2l3l4l5l6[0]*v_list[0]+l1l2l3l4l5l6[1]*v_list[1]+l1l2l3l4l5l6[2]*v_list[2]+l1l2l3l4l5l6[3]*v_list[3]+l1l2l3l4l5l6[4]*v_list[4]+l1l2l3l4l5l6[5]*v_list[5]; +} +/*}}}*/ +/*FUNCTION GetJacobianInvert {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetJacobianInvert" +void Penta::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3l4){ + + double Jdet; + const int NDOF3=3; + + /*Call Jacobian routine to get the jacobian:*/ + GetJacobian(Jinv, xyz_list, gauss_l1l2l3l4); + + /*Invert Jacobian matrix: */ + MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet); +} +/*}}}*/ +/*FUNCTION GetMatPar {{{1*/ +void* Penta::GetMatPar(){ + return matpar; +} +/*}}}*/ +/*FUNCTION GetNodalFunctions {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetNodalFunctions" +void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){ + + /*This routine returns the values of the nodal functions at the gaussian point.*/ + + l1l2l3l4l5l6[0]=gauss_l1l2l3l4[0]*(1-gauss_l1l2l3l4[3])/2.0; + + l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0; + + l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0; + + l1l2l3l4l5l6[3]=gauss_l1l2l3l4[0]*(1+gauss_l1l2l3l4[3])/2.0; + + l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0; + + l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0; + +} +/*}}}*/ +/*FUNCTION GetNodes {{{1*/ +void Penta::GetNodes(void** vpnodes){ int i; - const int calculationdof=3; - const int numgrids=6; - int DOFPERGRID=1; - - /*Same thing in the basic coordinate system: */ - double dh1dh6_basic[calculationdof][numgrids]; - - /*Get dh1dh2dh3 in basic coordinates system : */ - GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord); - - /*Build B': */ - for (i=0;iCreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type); - delete tria; - return; - } -} -/*}}}*/ -/*FUNCTION Penta CreatePVectorThermal {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::CreatePVectorThermal" -void Penta::CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ - - - /*indexing: */ - int i,j; - int found=0; - - const int numgrids=6; - const int NDOF1=1; - const int numdof=numgrids*NDOF1; - int doflist[numdof]; - int numberofdofspernode; - - /*Grid data: */ - double xyz_list[numgrids][3]; - - /* gaussian points: */ - int num_area_gauss,igarea,igvert; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* vert_gauss_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - double gauss_weight,area_gauss_weight,vert_gauss_weight; - double gauss_coord[4]; - int area_order=2; - int num_vert_gauss=3; - - double dt; - double vx_list[numgrids]; - double vy_list[numgrids]; - double vz_list[numgrids]; - double vxvyvz_list[numgrids][3]; - double temperature_list[numgrids]; - double temperature; - - /*Material properties: */ - double gravity,rho_ice,rho_water; - double mixed_layer_capacity,heatcapacity; - double beta,meltingpoint,thermal_exchange_velocity; - - /* element parameters: */ - int friction_type; - - int dofs[3]={0,1,2}; - int dofs1[1]={0}; - - /*matrices: */ - double P_terms[numdof]={0.0}; - double L[numdof]; - double l1l2l3[3]; - double alpha2_list[3]; - double basalfriction_list[3]={0.0}; - double basalfriction; - double epsilon[6]; - double epsilon_sqr[3][3]; - double epsilon_matrix[3][3]; - - double Jdet; - double viscosity; - double epsilon_eff; - double phi; - double t_pmp; - double scalar; - double scalar_def; - double scalar_ocean; - double scalar_transient; - - /*Collapsed formulation: */ - Tria* tria=NULL; - ParameterInputs* inputs=NULL; - - /*If on water, skip: */ - if(onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*recovre material parameters: */ - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - heatcapacity=matpar->GetHeatCapacity(); - beta=matpar->GetBeta(); - meltingpoint=matpar->GetMeltingPoint(); - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - if(dt){ - found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); - } - - for(i=0;iGetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /* Get nodal functions */ - GetNodalFunctions(&L[0], gauss_coord); - - /*Build deformational heating: */ - GetPhi(&phi, &epsilon[0], viscosity); - - /*Build pe_gaussian */ - scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; - if(dt){ - scalar_def=scalar_def*dt; - } - - for(i=0;iCreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type); - delete tria; - } - - /* Geothermal flux on ice sheet base and basal friction */ - if(onbed && !shelf){ - - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type); - delete tria; - } - extern int my_rank; - -cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&vert_gauss_coord); - xfree((void**)&area_gauss_weights); - xfree((void**)&vert_gauss_weights); - -} -/*}}}*/ -/*FUNCTION Penta CreatePVectorMelting {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Penta::CreatePVectorMelting" -void Penta::CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ - return; -} -/*}}}*/ -/*FUNCTION Penta GetPhi {{{1*/ + Node** pnodes=(Node**)vpnodes; + + for(i=0;i<6;i++){ + pnodes[i]=nodes[i]; + } +} +/*}}}*/ +/*FUNCTION GetOnBed {{{1*/ +int Penta::GetOnBed(){ + return onbed; +} +/*}}}*/ +/*FUNCTION GetPhi {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::GetPhi" @@ -4095,5 +3616,199 @@ } /*}}}*/ -/*FUNCTION Penta MassFlux {{{1*/ +/*FUNCTION GetShelf {{{1*/ +int Penta::GetShelf(){ + return shelf; +} +/*}}}*/ +/*FUNCTION GetStrainRate {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetStrainRate" +void Penta::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4){ + + int i; + const int numgrids=6; + const int NDOF2=2; + + double B[5][NDOF2*numgrids]; + + /*Get B matrix: */ + GetB(&B[0][0], xyz_list, gauss_l1l2l3l4); + +#ifdef _ISSM_DEBUG_ + printf("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); + printf(" [ %lf %lf \n",B[1][0],B[1][1]); + printf(" [ %lf %lf ]\n",B[2][0],B[2][1]); + printf(" [ %lf %lf ]\n",B[3][0],B[3][1]); + printf(" [ %lf %lf ]\n",B[4][0],B[4][1]); + + printf("B for grid2 : [ %lf %lf \n",B[0][2],B[0][3]); + printf(" [ %lf %lf \n",B[1][2],B[1][3]); + printf(" [ %lf %lf ]\n",B[2][2],B[2][3]); + printf(" [ %lf %lf ]\n",B[3][2],B[3][3]); + printf(" [ %lf %lf ]\n",B[4][2],B[4][3]); + + printf("B for grid3 : [ %lf %lf \n", B[0][4],B[0][5]); + printf(" [ %lf %lf \n", B[1][4],B[1][5]); + printf(" [ %lf %lf ]\n",B[2][4],B[2][5]); + printf(" [ %lf %lf ]\n",B[3][4],B[3][5]); + printf(" [ %lf %lf ]\n",B[4][4],B[4][5]); + + printf("B for grid4 : [ %lf %lf \n", B[0][6],B[0][7]); + printf(" [ %lf %lf \n", B[1][6],B[1][7]); + printf(" [ %lf %lf ]\n",B[2][6],B[2][7]); + printf(" [ %lf %lf ]\n",B[3][6],B[3][7]); + printf(" [ %lf %lf ]\n",B[4][6],B[4][7]); + + printf("B for grid5 : [ %lf %lf \n", B[0][8],B[0][9]); + printf(" [ %lf %lf \n", B[1][8],B[1][9]); + printf(" [ %lf %lf ]\n",B[2][8],B[2][9]); + printf(" [ %lf %lf ]\n",B[3][8],B[3][9]); + printf(" [ %lf %lf ]\n",B[4][8],B[4][9]); + + printf("B for grid6 : [ %lf %lf \n", B[0][10],B[0][11]); + printf(" [ %lf %lf \n", B[1][10],B[1][11]); + printf(" [ %lf %lf ]\n",B[2][10],B[2][11]); + printf(" [ %lf %lf ]\n",B[3][10],B[3][11]); + printf(" [ %lf %lf ]\n",B[4][10],B[4][11]); + +#endif + + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, + velocity,NDOF2*numgrids,1,0, + epsilon,0); + +} +/*}}}*/ +/*FUNCTION GetStrainRateStokes {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetStrainRateStokes" +void Penta::GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord){ + + int i,j; + + const int numgrids=6; + const int DOFVELOCITY=3; + double B[8][27]; + double B_reduced[numgrids][DOFVELOCITY*numgrids]; + + /*Get B matrix: */ + GetBStokes(&B[0][0], xyz_list, gauss_coord); + + /*Create a reduced matrix of B to get rid of pressure */ + for (i=0;i<6;i++){ + for (j=0;j<3;j++){ + B_reduced[i][j]=B[i][j]; + } + for (j=4;j<7;j++){ + B_reduced[i][j-1]=B[i][j]; + } + for (j=8;j<11;j++){ + B_reduced[i][j-2]=B[i][j]; + } + for (j=12;j<15;j++){ + B_reduced[i][j-3]=B[i][j]; + } + for (j=16;j<19;j++){ + B_reduced[i][j-4]=B[i][j]; + } + for (j=20;j<23;j++){ + B_reduced[i][j-5]=B[i][j]; + } + } + /*Multiply B by velocity, to get strain rate: */ + MatrixMultiply( &B_reduced[0][0],6,DOFVELOCITY*numgrids, 0, velocity,DOFVELOCITY*numgrids,1,0,epsilon, 0); +} +/*}}}*/ +/*FUNCTION GetThicknessList {{{1*/ +void Penta::GetThicknessList(double* thickness_list){ + + int i; + for(i=0;i<6;i++)thickness_list[i]=h[i]; +} +/*}}}*/ +/*FUNCTION Gradj {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::Gradj" +void Penta::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ + + /*If on water, skip grad (=0): */ + if(onwater)return; + + if (strcmp(control_type,"drag")==0){ + GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type); + } + else if (strcmp(control_type,"B")==0){ + GradjB( grad_g, inputs,analysis_type,sub_analysis_type); + } + else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type)); +} +/*}}}*/ +/*FUNCTION GradjDrag {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GradjDrag" +void Penta::GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ + + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*If on shelf, skip: */ + if(shelf)return; + + /*Bail out if this element does not touch the bedrock: */ + if (!onbed) return; + + if (sub_analysis_type==HorizAnalysisEnum()){ + + /*MacAyeal or Pattyn*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type); + delete tria; + return; + } + else if (sub_analysis_type==StokesAnalysisEnum()){ + + /*Stokes*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type); + delete tria; + return; + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); +} +/*}}}*/ +/*FUNCTION GradjB {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::GradjB" +void Penta::GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ + + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + if (collapse){ + /*Bail out element if collapsed (2d) and not on bed*/ + if (!onbed) return; + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute gardj*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type); + delete tria; + return; + } + else{ + /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type); + delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION MassFlux {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::MassFlux" @@ -4102,2 +3817,292 @@ } /*}}}*/ +/*FUNCTION Misfit {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::Misfit" +double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){ + + double J; + Tria* tria=NULL; + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return 0; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute Misfit*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->Misfit(inputs,analysis_type,sub_analysis_type); + delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->Misfit(inputs,analysis_type,sub_analysis_type); + delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION MyRank {{{1*/ +int Penta::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ +/*FUNCTION ReduceMatrixStokes {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "ReduceMatrixStokes" +void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){ + + int i,j; + + double Kii[24][24]; + double Kib[24][3]; + double Kbb[3][3]; + double Kbi[3][24]; + double Kbbinv[3][3]; + double KibKbbinv[24][3]; + double Kright[24][24]; + + /*Create the four matrices used for reduction */ + for(i=0;i<24;i++){ + for(j=0;j<24;j++){ + Kii[i][j]=*(Ke_temp+27*i+j); + } + for(j=0;j<3;j++){ + Kib[i][j]=*(Ke_temp+27*i+j+24); + } + } + for(i=0;i<3;i++){ + for(j=0;j<24;j++){ + Kbi[i][j]=*(Ke_temp+27*(i+24)+j); + } + for(j=0;j<3;j++){ + Kbb[i][j]=*(Ke_temp+27*(i+24)+j+24); + } + } + + /*Inverse the matrix corresponding to bubble part Kbb */ + GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]); + + /*Multiply matrices to create the reduce matrix Ke_reduced */ + MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0); + MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Kbi[0][0],3,24,0,&Kright[0][0],0); + + /*Affect value to the reduced matrix */ + for(i=0;i<24;i++){ + for(j=0;j<24;j++){ + *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j]; + } + } +} +/*}}}*/ +/*FUNCTION ReduceVectorStokes {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::ReduceVectorStokes" +void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){ + + int i,j; + + double Pi[24]; + double Pb[3]; + double Kbb[3][3]; + double Kib[24][3]; + double Kbbinv[3][3]; + double KibKbbinv[24][3]; + double Pright[24]; + + /*Create the four matrices used for reduction */ + for(i=0;i<24;i++){ + Pi[i]=*(Pe_temp+i); + } + for(i=0;i<3;i++){ + Pb[i]=*(Pe_temp+i+24); + } + for(j=0;j<3;j++){ + for(i=0;i<24;i++){ + Kib[i][j]=*(Ke_temp+3*i+j); + } + for(i=0;i<3;i++){ + Kbb[i][j]=*(Ke_temp+3*(i+24)+j); + } + } + + /*Inverse the matrix corresponding to bubble part Kbb */ + GetMatrixInvert(&Kbbinv[0][0], &Kbb[0][0]); + + /*Multiply matrices to create the reduce matrix Ke_reduced */ + MatrixMultiply(&Kib[0][0],24,3,0,&Kbbinv[0][0],3,3,0,&KibKbbinv[0][0],0); + MatrixMultiply(&KibKbbinv[0][0],24,3,0,&Pb[0],3,1,0,&Pright[0],0); + + /*Affect value to the reduced matrix */ + for(i=0;i<24;i++){ + *(Pe_reduced+i)=Pi[i]-Pright[i]; + } +} +/*}}}*/ +/*FUNCTION SpawnTria {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::SpawnTria" +void* Penta::SpawnTria(int g0, int g1, int g2){ + + /*out of grids g0,g1 and g2 from Penta, build a tria element: */ + Tria* tria=NULL; + double tria_h[3]; + double tria_s[3]; + double tria_b[3]; + double tria_c[3]; + double tria_k[3]; + double tria_melting[3]; + double tria_accumulation[3]; + double tria_geothermalflux[3]; + + /*configuration: */ + int tria_node_ids[3]; + Node* tria_nodes[3]; + int tria_node_offsets[3]; + + tria_h[0]=h[g0]; + tria_h[1]=h[g1]; + tria_h[2]=h[g2]; + + tria_s[0]=s[g0]; + tria_s[1]=s[g1]; + tria_s[2]=s[g2]; + + tria_b[0]=b[g0]; + tria_b[1]=b[g1]; + tria_b[2]=b[g2]; + + tria_k[0]=k[g0]; + tria_k[1]=k[g1]; + tria_k[2]=k[g2]; + + tria_melting[0]=melting[g0]; + tria_melting[1]=melting[g1]; + tria_melting[2]=melting[g2]; + + tria_accumulation[0]=accumulation[g0]; + tria_accumulation[1]=accumulation[g1]; + tria_accumulation[2]=accumulation[g2]; + + tria_geothermalflux[0]=geothermalflux[g0]; + tria_geothermalflux[1]=geothermalflux[g1]; + tria_geothermalflux[2]=geothermalflux[g2]; + + tria_node_ids[0]=node_ids[g0]; + tria_node_ids[1]=node_ids[g1]; + tria_node_ids[2]=node_ids[g2]; + + tria_nodes[0]=nodes[g0]; + tria_nodes[1]=nodes[g1]; + tria_nodes[2]=nodes[g2]; + + tria_node_offsets[0]=node_offsets[g0]; + tria_node_offsets[1]=node_offsets[g1]; + tria_node_offsets[2]=node_offsets[g2]; + + tria= new Tria(id,mid,mparid,numparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,onwater); + + tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets); + tria->MaticeConfiguration(matice,matice_offset); + tria->MatparConfiguration(matpar,matpar_offset); + tria->NumparConfiguration(numpar,numpar_offset); + + return tria; + +} +/*}}}*/ +/*FUNCTION UpdateFromInputs {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::UpdateFromInputs" +void Penta::UpdateFromInputs(void* vinputs){ + + int dofs[1]={0}; + double temperature_list[6]; + double temperature_average; + double B_list[6]; + double B_average; + + ParameterInputs* inputs=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*Update internal data if inputs holds new values: */ + inputs->Recover("thickness",&h[0],1,dofs,6,(void**)nodes); + inputs->Recover("surface",&s[0],1,dofs,6,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,6,(void**)nodes); + inputs->Recover("drag",&k[0],1,dofs,6,(void**)nodes); + inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes); + inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes); + + //Update material if necessary + if(inputs->Recover("temperature",&temperature_list[0],1,dofs,6,(void**)nodes)){ + if(matice && !collapse){ + //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2]) + // +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0; + temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; + B_average=Paterson(temperature_average); + matice->SetB(B_average); + } + } + + if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){ + if(matice && collapse){ + temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; + B_average=Paterson(temperature_average); + //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2]) + // +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0; + matice->SetB(B_average); + } + } + + if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){ + if(matice){ + B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0; + matice->SetB(B_average); + } + } + +} +/*}}}*/ +/*FUNCTION SurfaceNormal {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::SurfaceNormal" +void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ + + int i; + double v13[3]; + double v23[3]; + double normal[3]; + double normal_norm; + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0][i]-xyz_list[2][i]; + v23[i]=xyz_list[1][i]-xyz_list[2][i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + + normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) ); + + *(surface_normal)=normal[0]/normal_norm; + *(surface_normal+1)=normal[1]/normal_norm; + *(surface_normal+2)=normal[2]/normal_norm; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 2712) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 2713) @@ -24,4 +24,5 @@ //#define _DEBUGGAUSS_ +/*Object constructors and destructor*/ /*FUNCTION Tria constructor {{{1*/ Tria::Tria(){ @@ -72,9 +73,2020 @@ } /*}}}*/ -/*FUNCTION Echo {{{1*/ + +/*Object marshall*/ +/*FUNCTION Marshall {{{1*/ +void Tria::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of Tria: */ + enum_type=TriaEnum(); + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall Tria data: */ + memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid); + memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid); + memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids); + memcpy(marshalled_dataset,&nodes,sizeof(nodes));marshalled_dataset+=sizeof(nodes); + memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets); + memcpy(marshalled_dataset,&matice,sizeof(matice));marshalled_dataset+=sizeof(matice); + memcpy(marshalled_dataset,&matice_offset,sizeof(matice_offset));marshalled_dataset+=sizeof(matice_offset); + memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar); + memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset); + memcpy(marshalled_dataset,&numparid,sizeof(numparid));marshalled_dataset+=sizeof(numparid); + memcpy(marshalled_dataset,&numpar,sizeof(numpar));marshalled_dataset+=sizeof(numpar); + memcpy(marshalled_dataset,&numpar_offset,sizeof(numpar_offset));marshalled_dataset+=sizeof(numpar_offset); + memcpy(marshalled_dataset,&h,sizeof(h));marshalled_dataset+=sizeof(h); + memcpy(marshalled_dataset,&s,sizeof(s));marshalled_dataset+=sizeof(s); + memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b); + memcpy(marshalled_dataset,&k,sizeof(k));marshalled_dataset+=sizeof(k); + memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting); + memcpy(marshalled_dataset,&accumulation,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation); + memcpy(marshalled_dataset,&geothermalflux,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux); + memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type); + memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed); + memcpy(marshalled_dataset,&onwater,sizeof(onwater));marshalled_dataset+=sizeof(onwater); + memcpy(marshalled_dataset,&p,sizeof(p));marshalled_dataset+=sizeof(p); + memcpy(marshalled_dataset,&q,sizeof(q));marshalled_dataset+=sizeof(q); + memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf); + + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION MarshallSize {{{1*/ +int Tria::MarshallSize(){ + return sizeof(id) + +sizeof(mid) + +sizeof(mparid) + +sizeof(node_ids) + +sizeof(nodes) + +sizeof(node_offsets) + +sizeof(matice) + +sizeof(matice_offset) + +sizeof(matpar) + +sizeof(matpar_offset) + +sizeof(numparid) + +sizeof(numpar) + +sizeof(numpar_offset) + +sizeof(h) + +sizeof(s) + +sizeof(b) + +sizeof(k) + +sizeof(melting) + +sizeof(accumulation) + +sizeof(geothermalflux) + +sizeof(friction_type) + +sizeof(onbed) + +sizeof(onwater) + +sizeof(p) + +sizeof(q) + +sizeof(shelf) + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION Demarshall {{{1*/ +void Tria::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + + memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid); + memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid); + memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids); + memcpy(&nodes,marshalled_dataset,sizeof(nodes));marshalled_dataset+=sizeof(nodes); + memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets); + memcpy(&matice,marshalled_dataset,sizeof(matice));marshalled_dataset+=sizeof(matice); + memcpy(&matice_offset,marshalled_dataset,sizeof(matice_offset));marshalled_dataset+=sizeof(matice_offset); + memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar); + memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset); + memcpy(&numparid,marshalled_dataset,sizeof(numparid));marshalled_dataset+=sizeof(numparid); + memcpy(&numpar,marshalled_dataset,sizeof(numpar));marshalled_dataset+=sizeof(numpar); + memcpy(&numpar_offset,marshalled_dataset,sizeof(numpar_offset));marshalled_dataset+=sizeof(numpar_offset); + memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h); + memcpy(&s,marshalled_dataset,sizeof(s));marshalled_dataset+=sizeof(s); + memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b); + memcpy(&k,marshalled_dataset,sizeof(k));marshalled_dataset+=sizeof(k); + memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting); + memcpy(&accumulation,marshalled_dataset,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation); + memcpy(&geothermalflux,marshalled_dataset,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux); + memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type); + memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed); + memcpy(&onwater,marshalled_dataset,sizeof(onwater));marshalled_dataset+=sizeof(onwater); + memcpy(&p,marshalled_dataset,sizeof(p));marshalled_dataset+=sizeof(p); + memcpy(&q,marshalled_dataset,sizeof(q));marshalled_dataset+=sizeof(q); + memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf); + + /*nodes and materials are not pointing to correct objects anymore:*/ + for(i=0;i<3;i++)nodes[i]=NULL; + matice=NULL; + matpar=NULL; + numpar=NULL; + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ + +/*Object functions*/ +/*FUNCTION ComputePressure {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::ComputePressure" +void Tria::ComputePressure(Vec pg){ + + int i; + const int numgrids=3; + int doflist[numgrids]; + double pressure[numgrids]; + double rho_ice,g; + + /*Get dof list on which we will plug the pressure values: */ + GetDofList1(&doflist[0]); + + /*pressure is lithostatic: */ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + for(i=0;iRecover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iGetRhoIce()); + printf(" gravity: %g \n",matpar->GetG()) + printf(" rho_water: %g \n",matpar->GetRhoWater()); + printf(" Velocity: \n"); + for (i=0;iGetB()); + printf(" drag [%g %g %g ]\n",k[0],k[1],k[2]); + printf(" thickness [%g %g %g]\n",h[0],h[1],h[2]); + printf(" surface [%g %g %g ]\n",s[0],s[1],s[2]); + printf(" bed [%g %g %g]\n",b[0],b[1],b[2]); + } +#endif + + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + +#ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;iGetViscosity2d(&viscosity, &epsilon[0]); + matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=newviscosity*thickness*gauss_weight*Jdet; + + for (i=0;i<3;i++){ + D[i][i]=D_scalar; + } + +#ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian loop %i\n",ig); + printf(" thickness %g\n",thickness); + printf(" slope [%g,%g]\n",slope[0],slope[1]); + printf(" alpha2_list [%g,%g,%g]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); + printf(" epsilon [%g,%g,%g]\n",epsilon[0],epsilon[1],epsilon[2]); + printf(" Matice: \n"); + matice->Echo(); + printf(" Matpar: \n"); + matpar->Echo(); + printf("\n viscosity: %g \n",viscosity); + printf(" jacobian: %g \n",Jdet); + printf(" gauss_weight: %g \n",gauss_weight); + } +#endif + + /*Get B and Bprime matrices: */ + GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],3,numdof,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iRecover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + + /*Build alpha2_list used by drag stiffness matrix*/ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=matpar->GetG(); + friction->rho_ice=matpar->GetRhoIce(); + friction->rho_water=matpar->GetRhoWater(); + friction->K=&k[0]; + friction->bed=&b[0]; + friction->thickness=&h[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=p; + friction->q=q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); + } + #endif + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, + //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ + GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); + slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); + + if (slope_magnitude>MAXSLOPE){ + alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); + alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); + alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); + } + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get L matrix: */ + GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + + /*Now, take care of the basal friction if there is any: */ + GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); + + DL_scalar=alpha2*gauss_weight*Jdet; + for (i=0;i<2;i++){ + DL[i][i]=DL_scalar; + } + + /* Do the triple producte tL*D*L: */ + TripleMultiply( &L[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &L[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + for( i=0; iGetLatentHeat(); + heatcapacity=matpar->GetHeatCapacity(); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights: */ + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igRecover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity_average in inputs!"); + + for(i=0;iRecover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //Create Artificial diffusivity once for all if requested + if(numpar->artdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); + v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); + + K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); + K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igartdiff){ + + /* Compute artificial diffusivity */ + KDL[0][0]=DL_scalar*K[0][0]; + KDL[1][1]=DL_scalar*K[1][1]; + + TripleMultiply( &Bprime[0][0],2,numdof,1, + &KDL[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add artificial diffusivity matrix */ + for( i=0; iRecover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //recover material parameters + mixed_layer_capacity=matpar->GetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + + + GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igRecover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity"); + inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set pe_g to 0: */ + for(i=0;iGetG()); + printf("rho_ice %g\n",matpar->GetRhoIce()); + printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]); + printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]); + printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]); + printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]); + } + #endif + + + /* Get gaussian points and weights: */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;iGetRhoIce()*matpar->GetG()*thickness; + + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian %i\n",ig); + printf(" thickness %g\n",thickness); + printf(" slope(%g,%g)\n",slope[0],slope[1]); + printf(" Jdet %g\n",Jdet); + printf(" gaussweigth %g\n",gauss_weight); + printf(" l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]); + if(friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); + } + #endif + + /*Build pe_g_gaussian vector: */ + if(friction_type==1){ + for (i=0;iterms\n",ig); + for( i=0; inrows; i++){ + printf("%g ",*(pe_g->terms+i)); + } + printf("\n"); + printf(" pe_g->row_indices\n",ig); + for( i=0; inrows; i++){ + printf("%i ",*(pe_g->row_indices+i)); + } + } + #endif + + /*Add pe_g to global vector pg: */ + VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + +} +/*}}}*/ +/*FUNCTION CreatePVectorPrognostic {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::CreatePVectorPrognostic" +void Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ + + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* matrix */ + double pe_g[numgrids]={0.0}; + double L[numgrids]; + double Jdettria; + + /*input parameters for structural analysis (diagnostic): */ + double accumulation_list[numgrids]={0.0}; + double accumulation_g; + double melting_list[numgrids]={0.0}; + double melting_g; + double thickness_list[numgrids]={0.0}; + double thickness_g; + double dt; + int dofs[1]={0}; + int found=0; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover extra inputs from users, at current convergence iteration: */ + found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!"); + found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!"); + found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find thickness in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + beta=matpar->GetBeta(); + meltingpoint=matpar->GetMeltingPoint(); + + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); + + /* Ice/ocean heat exchange flux on ice shelf base */ + + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + for (ig=0; igGetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + for(i=0;ielement_type=(char*)xmalloc((strlen("3d")+1)*sizeof(char)); + strcpy(friction->element_type,"3d"); + + friction->gravity=matpar->GetG(); + friction->rho_ice=matpar->GetRhoIce(); + friction->rho_water=matpar->GetRhoWater(); + friction->K=&k[0]; + friction->bed=&b[0]; + friction->thickness=&h[0]; + friction->velocities=&vxvyvz_list[0][0]; + friction->p=p; + friction->q=q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + /* Compute basal friction */ + for(i=0;iRecover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); + if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"missing velocity_obs input parameter"); + } + if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"missing velocity input parameter"); + } + + for(i=0;imeanvel/(obs_vx_list[i]+numpar->epsvel),2); + scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),2); + if(obs_vx_list[i]==0)scalex=0; + if(obs_vy_list[i]==0)scaley=0; + relativex_list[i]=scalex*(obs_vx_list[i]-vx_list[i]); + relativey_list[i]=scaley*(obs_vy_list[i]-vy_list[i]); + } + } + else if(fit==2){ + /*We are using a logarithmic misfit: */ + for (i=0;iepsvel; //epsvel to avoid velocity being nil. + obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. + scale=-8*pow(numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + logarithmicx_list[i]=scale*vx_list[i]; + logarithmicy_list[i]=scale*vy_list[i]; + } + } + else{ + /*Not supported yet! : */ + throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit)); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + +#ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;iGetArea(); + + /*Get xyz list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + x1=xyz_list[0][0]; y1=xyz_list[0][1]; + x2=xyz_list[1][0]; y2=xyz_list[1][1]; + x3=xyz_list[2][0]; y3=xyz_list[2][1]; + + if(which_one==1){ + /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ + return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; + } + else if(which_one==2){ + /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ + return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; + } + else if(which_one==3){ + /*Get third area coordinate 1-area1-area2: */ + return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!")); +} +/*}}}*/ +/*FUNCTION GetB {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::GetB" + +void Tria::GetB(double* B, double* xyz_list, double* gauss_l1l2l3){ + + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + * For grid i, Bi can be expressed in the basic coordinate system + * by: + * Bi_basic=[ dh/dx 0 ] + * [ 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 3x(NDOF2*numgrids) + */ + int i; - - DataSet* loadsin=NULL; - DataSet* nodesin=NULL; - DataSet* materialsin=NULL; - DataSet* parametersin=NULL; - - /*Recover pointers :*/ - loadsin=(DataSet*)ploadsin; - nodesin=(DataSet*)pnodesin; - materialsin=(DataSet*)pmaterialsin; - parametersin=(DataSet*)pparametersin; - - /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */ - ResolvePointers((Object**)nodes,node_ids,node_offsets,3,nodesin); - - /*Same for materials: */ - ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin); - ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin); - - /*Same for numpar: */ - ResolvePointers((Object**)&numpar,&numparid,&numpar_offset,1,parametersin); - -} -/*}}}*/ -/*FUNCTION CreateKMatrix {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreateKMatrix" - -void Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==ControlAnalysisEnum()){ - - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); - } - else if (analysis_type==DiagnosticAnalysisEnum()){ - - if (sub_analysis_type==HorizAnalysisEnum()){ - - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); - } - else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); - - } - else if (analysis_type==SlopeComputeAnalysisEnum()){ - - CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type); - - } - else if (analysis_type==PrognosticAnalysisEnum()){ - - CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type); - - } - else{ - throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); - } - -} -/*}}}*/ -/*FUNCTION CreateKMatrixDiagnosticHoriz {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz" - -void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ - - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=3; - const int numdof=2*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - /* material data: */ - double viscosity; //viscosity - double newviscosity; //viscosity - double oldviscosity; //viscosity - - /* strain rate: */ - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - double oldepsilon[3]; /* oldepsilon=[exx,eyy,exy];*/ - - /* matrices: */ - double B[3][numdof]; - double Bprime[3][numdof]; - double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. - double D_scalar; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - - double Jdet; - - /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - double oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - double thickness; - int dofs[2]={0,1}; - - ParameterInputs* inputs=NULL; - - /*First, if we are on water, return empty matrix: */ - if(onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;iGetRhoIce()); - printf(" gravity: %g \n",matpar->GetG()) - printf(" rho_water: %g \n",matpar->GetRhoWater()); - printf(" Velocity: \n"); - for (i=0;iGetB()); - printf(" drag [%g %g %g ]\n",k[0],k[1],k[2]); - printf(" thickness [%g %g %g]\n",h[0],h[1],h[2]); - printf(" surface [%g %g %g ]\n",s[0],s[1],s[2]); - printf(" bed [%g %g %g]\n",b[0],b[1],b[2]); + const int NDOF2=2; + const int numgrids=3; + + double dh1dh2dh3_basic[NDOF2][numgrids]; + + + /*Get dh1dh2dh3 in basic coordinate system: */ + GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3); + + #ifdef _ISSM_DEBUG_ + for (i=0;i<3;i++){ + printf("Node %i dh/dx=%lf dh/dy=%lf \n",i,dh1dh2dh3_basic[0][i],dh1dh2dh3_basic[1][i]); } #endif - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;iGetViscosity2d(&viscosity, &epsilon[0]); - matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ - newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*thickness*gauss_weight*Jdet; - - for (i=0;i<3;i++){ - D[i][i]=D_scalar; - } - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian loop %i\n",ig); - printf(" thickness %g\n",thickness); - printf(" slope [%g,%g]\n",slope[0],slope[1]); - printf(" alpha2_list [%g,%g,%g]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); - printf(" epsilon [%g,%g,%g]\n",epsilon[0],epsilon[1],epsilon[2]); - printf(" Matice: \n"); - matice->Echo(); - printf(" Matpar: \n"); - matpar->Echo(); - printf("\n viscosity: %g \n",viscosity); - printf(" jacobian: %g \n",Jdet); - printf(" gauss_weight: %g \n",gauss_weight); - } - #endif - - /*Get B and Bprime matrices: */ - GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iRecover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find velocity_average in inputs!"); - - for(i=0;iRecover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - //Create Artificial diffusivity once for all if requested - if(numpar->artdiff){ - //Get the Jacobian determinant - gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); - - //Build K matrix (artificial diffusivity matrix) - v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); - v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); - - K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); - K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igartdiff){ - - /* Compute artificial diffusivity */ - KDL[0][0]=DL_scalar*K[0][0]; - KDL[1][1]=DL_scalar*K[1][1]; - - TripleMultiply( &Bprime[0][0],2,numdof,1, - &KDL[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add artificial diffusivity matrix */ - for( i=0; iRecover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - - /*Build alpha2_list used by drag stiffness matrix*/ - Friction* friction=NewFriction(); - - /*Initialize all fields: */ - friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); - strcpy(friction->element_type,"2d"); - - friction->gravity=matpar->GetG(); - friction->rho_ice=matpar->GetRhoIce(); - friction->rho_water=matpar->GetRhoWater(); - friction->K=&k[0]; - friction->bed=&b[0]; - friction->thickness=&h[0]; - friction->velocities=&vxvy_list[0][0]; - friction->p=p; - friction->q=q; - - /*Compute alpha2_list: */ - FrictionGetAlpha2(&alpha2_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); - } - #endif - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, - //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ - GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); - slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - - if (slope_magnitude>MAXSLOPE){ - alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); - } - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - - /*Now, take care of the basal friction if there is any: */ - GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); - - DL_scalar=alpha2*gauss_weight*Jdet; - for (i=0;i<2;i++){ - DL[i][i]=DL_scalar; - } - - /* Do the triple producte tL*D*L: */ - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; iGetG()); - printf("rho_ice %g\n",matpar->GetRhoIce()); - printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]); - printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]); - printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]); - printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]); - } - #endif - - - /* Get gaussian points and weights: */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;iGetRhoIce()*matpar->GetG()*thickness; - - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian %i\n",ig); - printf(" thickness %g\n",thickness); - printf(" slope(%g,%g)\n",slope[0],slope[1]); - printf(" Jdet %g\n",Jdet); - printf(" gaussweigth %g\n",gauss_weight); - printf(" l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]); - if(friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); - } - #endif - - /*Build pe_g_gaussian vector: */ - if(friction_type==1){ - for (i=0;iterms\n",ig); - for( i=0; inrows; i++){ - printf("%g ",*(pe_g->terms+i)); - } - printf("\n"); - printf(" pe_g->row_indices\n",ig); - for( i=0; inrows; i++){ - printf("%i ",*(pe_g->row_indices+i)); - } - } - #endif - - /*Add pe_g to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - - cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - -} -/*}}}*/ -/*FUNCTION CreatePVectorPrognostic {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreatePVectorPrognostic" -void Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ - - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - /* matrix */ - double pe_g[numgrids]={0.0}; - double L[numgrids]; - double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ - double accumulation_list[numgrids]={0.0}; - double accumulation_g; - double melting_list[numgrids]={0.0}; - double melting_g; - double thickness_list[numgrids]={0.0}; - double thickness_g; - double dt; - int dofs[1]={0}; - int found=0; - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!"); - found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!"); - found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find thickness in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igRecover("thickness",&h[0],1,dofs,3,(void**)nodes); - inputs->Recover("surface",&s[0],1,dofs,3,(void**)nodes); - inputs->Recover("bed",&b[0],1,dofs,3,(void**)nodes); - inputs->Recover("drag",&k[0],1,dofs,3,(void**)nodes); - inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes); - inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes); - inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes); - - //Update material if necessary - if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){ - temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2])/3.0; - B_average=Paterson(temperature_average); - matice->SetB(B_average); - } - - if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){ - B_average=(B_list[0]+B_list[1]+B_list[2])/3.0; - matice->SetB(B_average); - } - + /*Build B: */ + for (i=0;iRecover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); - if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"missing velocity_obs input parameter"); - } + GetDofList1(&doflist1[0]); + + /* Set grade_g to 0: */ + for(i=0;iRecover("thickness",&h[0],1,dofs,numgrids,(void**)nodes); if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ throw ErrorException(__FUNCT__,"missing velocity input parameter"); } - + if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"missing adjoint input parameter"); + } + if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"parameter B not found in input"); + } + + /*Initialize parameter lists: */ for(i=0;imeanvel/(obs_vx_list[i]+numpar->epsvel),2); - scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),2); - if(obs_vx_list[i]==0)scalex=0; - if(obs_vy_list[i]==0)scaley=0; - relativex_list[i]=scalex*(obs_vx_list[i]-vx_list[i]); - relativey_list[i]=scaley*(obs_vy_list[i]-vy_list[i]); - } - } - else if(fit==2){ - /*We are using a logarithmic misfit: */ - for (i=0;iepsvel; //epsvel to avoid velocity being nil. - obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. - scale=-8*pow(numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); - logarithmicx_list[i]=scale*vx_list[i]; - logarithmicy_list[i]=scale*vy_list[i]; - } - } - else{ - /*Not supported yet! : */ - throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit)); - } - + adjx_list[i]=adjxadjy_list[i][0]; + adjy_list[i]=adjxadjy_list[i][1]; + } + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;iGetViscosityComplement(&viscosity_complement, &epsilon[0]); + + /*Get dvx, dvy, dadjx and dadjx: */ + GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); + GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); + GetParameterDerivativeValue(&dadjx[0], &adjx_list[0],&xyz_list[0][0], gauss_l1l2l3); + GetParameterDerivativeValue(&dadjy[0], &adjy_list[0],&xyz_list[0][0], gauss_l1l2l3); + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - #ifdef _ISSM_DEBUG_ - printf("Element id %i Jacobian determinant: %g\n",GetId(),Jdet); - #endif - + /* Get nodal functions value at gaussian point:*/ GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */ - if(fit==0){ - /*We are using an absolute misfit: */ - - /*Compute absolute(x/y) at gaussian point: */ - GetParameterValue(&absolutex, &absolutex_list[0],gauss_l1l2l3); - GetParameterValue(&absolutey, &absolutey_list[0],gauss_l1l2l3); - - /*compute Du*/ - for (i=0;icm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dB[0]+dh1dh2dh3_basic[1][i]*dB[1]); + + //min dampening + if(B_gausscm_mindmp_value){ + grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; } - } - else if(fit==1){ - /*We are using a relative misfit: */ - - /*Compute relative(x/y) at gaussian point: */ - GetParameterValue(&relativex, &relativex_list[0],gauss_l1l2l3); - GetParameterValue(&relativey, &relativey_list[0],gauss_l1l2l3); - - /*compute Du*/ - for (i=0;inumpar->cm_maxdmp_value){ + grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; } - } - else if(fit==2){ - /*We are using a logarithmic misfit: */ - - /*Compute logarithmic(x/y) at gaussian point: */ - GetParameterValue(&logarithmicx, &logarithmicx_list[0],gauss_l1l2l3); - GetParameterValue(&logarithmicy, &logarithmicy_list[0],gauss_l1l2l3); - - /*compute Du*/ - for (i=0;iRecover("thickness",&h[0],1,dofs,numgrids,(void**)nodes); - if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"missing velocity input parameter"); - } - if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"missing adjoint input parameter"); - } - if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"parameter B not found in input"); - } - - /*Initialize parameter lists: */ - for(i=0;iGetViscosityComplement(&viscosity_complement, &epsilon[0]); - - /*Get dvx, dvy, dadjx and dadjx: */ - GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dadjx[0], &adjx_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dadjy[0], &adjy_list[0],&xyz_list[0][0], gauss_l1l2l3); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - #ifdef _ISSM_DEBUG_ - printf("viscositycomp %g thickness %g dvx [%g %g] dvy [%g %g] dadjx [%g %g] dadjy[%g %g]\n",viscosity_complement,thickness,dvx[0],dvx[1],dvy[0],dvy[1],dadjx[0],dadjx[1],dadjy[0],dadjy[1]); - #endif - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],&xyz_list[0][0],gauss_l1l2l3); - - /*Get B derivative: dB/dx */ - GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3); - - /*Build gradje_g_gaussian vector (actually -dJ/dB): */ - for (i=0;icm_noisedmp*Jdet*gauss_weight*(dh1dh2dh3_basic[0][i]*dB[0]+dh1dh2dh3_basic[1][i]*dB[1]); - - //min dampening - if(B_gausscm_mindmp_value){ - grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; - } - - //max dampening - if(B_gauss>numpar->cm_maxdmp_value){ - grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; - } - - } - - /*Add grade_g_gaussian to grade_g: */ - for( i=0; iRecover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); - if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"missing velocity_obs input parameter"); - } - if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"missing velocity input parameter"); - } - - /*Initialize velocities: */ - for(i=0;imeanvel/(obs_vx_list[i]+numpar->epsvel),(double)2); - scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),(double)2); - if(obs_vx_list[i]==0)scalex=0; - if(obs_vy_list[i]==0)scaley=0; - relative_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2)); - } - } - else if(fit==2){ - /*We are using a logarithmic misfit: */ - for (i=0;iepsvel; //epsvel to avoid velocity being nil. - obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. - logarithmic_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); - } - } - else{ - /*Not supported yet! : */ - throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit)); - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - #ifdef _ISSM_DEBUG_ - for (i=0;icontrol_type,"drag")==0){ - if (!shelf){ - - //noise dampening - GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); - Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; - - } - } - else if (strcmp(numpar->control_type,"B")==0){ - if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ - throw ErrorException(__FUNCT__,"parameter B not found in input"); - } - //noise dampening - GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); - Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; - - //min dampening - GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3); - if(B_gausscm_mindmp_value){ - Jelem+=numpar->cm_mindmp_slope*B_gauss*Jdet*gauss_weight; - } - - //max dampening - if(B_gauss>numpar->cm_maxdmp_value){ - Jelem+=numpar->cm_maxdmp_slope*B_gauss*Jdet*gauss_weight; - } - } - else{ - throw ErrorException(__FUNCT__,exprintf("%s%s","unsupported control type: ",numpar->control_type)); - } - - /*Differents misfits are allowed: */ - if(fit==0){ - /*Compute absolute misfit at gaussian point: */ - GetParameterValue(&absolute, &absolute_list[0],gauss_l1l2l3); - - /*compute Misfit*/ - Jelem+=absolute*Jdet*gauss_weight; - } - else if(fit==1){ - /*Compute relative misfit at gaussian point: */ - GetParameterValue(&relative, &relative_list[0],gauss_l1l2l3); - - /*compute Misfit*/ - Jelem+=relative*Jdet*gauss_weight; - } - else if(fit==2){ - /*Compute logarithmic misfit at gaussian point: */ - GetParameterValue(&logarithmic, &logarithmic_list[0],gauss_l1l2l3); - - /*compute Misfit*/ - Jelem+=logarithmic*Jdet*gauss_weight; - } - else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!")); - - } - cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - - /*Return: */ - return Jelem; -} -/*}}}*/ -/*FUNCTION NodeConfiguration {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::NodeConfiguration" -void Tria::NodeConfiguration(int* tria_node_ids,Node* tria_nodes[3],int* tria_node_offsets){ - - int i; - for(i=0;i<3;i++){ - node_ids[i]=tria_node_ids[i]; - nodes[i]=tria_nodes[i]; - node_offsets[i]=tria_node_offsets[i]; - } - -} -/*}}}*/ -/*FUNCTION MaticeConfiguration {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::MaticeConfiguration" -void Tria::MaticeConfiguration(Matice* tria_matice,int tria_matice_offset){ - matice=tria_matice; - matice_offset=tria_matice_offset; -} -/*}}}*/ -/*FUNCTION MatparConfiguration {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::MatparConfiguration" -void Tria::MatparConfiguration(Matpar* tria_matpar,int tria_matpar_offset){ - - matpar=tria_matpar; - matpar_offset=tria_matpar_offset; - -} -/*}}}*/ -/*FUNCTION NumparConfiguration {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::NumparConfiguration" -void Tria::NumparConfiguration(Numpar* tria_numpar,int tria_numpar_offset){ - - numpar=tria_numpar; - numpar_offset=tria_numpar_offset; - -} -/*}}}*/ -/*FUNCTION CreateKMatrixDiagnosticSurfaceVert {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert" -void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ - - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - - /* surface normal: */ - double x4,y4,z4; - double x5,y5,z5; - double x6,y6,z6; - double v46[3]; - double v56[3]; - double normal[3]; - double norm_normal; - double nz; - - /*Matrices: */ - double DL_scalar; - double L[3]; - double Jdet; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;iRecover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity"); - inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set pe_g to 0: */ - for(i=0;iGetRhoIce(); - g=matpar->GetG(); - for(i=0;iRecover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - //recover material parameters - mixed_layer_capacity=matpar->GetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - - - GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igGetLatentHeat(); - heatcapacity=matpar->GetHeatCapacity(); - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Get gaussian points and weights: */ - GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igGetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - beta=matpar->GetBeta(); - meltingpoint=matpar->GetMeltingPoint(); - - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); - - /* Ice/ocean heat exchange flux on ice shelf base */ - - GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - for (ig=0; igGetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("dt",&dt); - if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); - - for(i=0;ielement_type=(char*)xmalloc((strlen("3d")+1)*sizeof(char)); - strcpy(friction->element_type,"3d"); - - friction->gravity=matpar->GetG(); - friction->rho_ice=matpar->GetRhoIce(); - friction->rho_water=matpar->GetRhoWater(); - friction->K=&k[0]; - friction->bed=&b[0]; - friction->thickness=&h[0]; - friction->velocities=&vxvyvz_list[0][0]; - friction->p=p; - friction->q=q; - - /*Compute alpha2_list: */ - FrictionGetAlpha2(&alpha2_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); - - /* Compute basal friction */ - for(i=0;imatpar->GetRhoIce(); @@ -3906,5 +3644,5 @@ /*Recover segment node locations: */ x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); - + /*Get xyz list: */ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); @@ -3940,61 +3678,328 @@ mass_flux= rho_ice*length*( - (1.0/3.0*(h1-h2)*(vx1-vx2)+1.0/2.0*h2*(vx1-vx2)+1.0/2.0*(h1-h2)*vx2+h2*vx2)*normal[0]+ - (1.0/3.0*(h1-h2)*(vy1-vy2)+1.0/2.0*h2*(vy1-vy2)+1.0/2.0*(h1-h2)*vy2+h2*vy2)*normal[1] + (1.0/3.0*(h1-h2)*(vx1-vx2)+1.0/2.0*h2*(vx1-vx2)+1.0/2.0*(h1-h2)*vx2+h2*vx2)*normal[0]+ + (1.0/3.0*(h1-h2)*(vy1-vy2)+1.0/2.0*h2*(vy1-vy2)+1.0/2.0*(h1-h2)*vy2+h2*vy2)*normal[1] ); return mass_flux; } /*}}}*/ -/*FUNCTION GetArea {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::GetArea" -double Tria::GetArea(void){ - - double area=0; +/*FUNCTION MaticeConfiguration {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::MaticeConfiguration" +void Tria::MaticeConfiguration(Matice* tria_matice,int tria_matice_offset){ + matice=tria_matice; + matice_offset=tria_matice_offset; +} +/*}}}*/ +/*FUNCTION MatparConfiguration {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::MatparConfiguration" +void Tria::MatparConfiguration(Matpar* tria_matpar,int tria_matpar_offset){ + + matpar=tria_matpar; + matpar_offset=tria_matpar_offset; + +} +/*}}}*/ +/*FUNCTION Misfit {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::Misfit" +double Tria::Misfit(void* vinputs,int analysis_type,int sub_analysis_type){ + + int i; + + /* output: */ + double Jelem=0; + + /* node data: */ const int numgrids=3; - double xyz_list[numgrids][3]; - double x1,y1,x2,y2,x3,y3; - - /*Get xyz list: */ + const int numdof=2*numgrids; + const int NDOF2=2; + int dofs1[1]={0}; + int dofs2[2]={0,1}; + double xyz_list[numgrids][3]; + + /* grid data: */ + double vxvy_list[numgrids][2]; + double vx_list[numgrids]; + double vy_list[numgrids]; + double obs_vxvy_list[numgrids][2]; + double obs_vx_list[numgrids]; + double obs_vy_list[numgrids]; + double absolute_list[numgrids]; + double relative_list[numgrids]; + double logarithmic_list[numgrids]; + double B[numgrids]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + double k_gauss; + double B_gauss; + + /* parameters: */ + double velocity_mag,obs_velocity_mag; + double absolute,relative,logarithmic; + double dk[NDOF2]; + double dB[NDOF2]; + + /* Jacobian: */ + double Jdet; + + /*relative and logarithmic control method :*/ + double scalex=1; + double scaley=1; + double fit=-1; + + ParameterInputs* inputs=NULL; + + /*If on water, return 0: */ + if(onwater)return 0; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - x1=xyz_list[0][0]; y1=xyz_list[0][1]; - x2=xyz_list[1][0]; y2=xyz_list[1][1]; - x3=xyz_list[2][0]; y3=xyz_list[2][1]; - - return x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1; -} -/*}}}*/ -/*FUNCTION GetAreaCoordinate {{{1*/ -#undef __FUNCT__ -#define __FUNCT__ "Tria::GetAreaCoordinate" -double Tria::GetAreaCoordinate(double x, double y, int which_one){ - - double area=0; - const int numgrids=3; - double xyz_list[numgrids][3]; - double x1,y1,x2,y2,x3,y3; - - /*Get area: */ - area=this->GetArea(); - - /*Get xyz list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - x1=xyz_list[0][0]; y1=xyz_list[0][1]; - x2=xyz_list[1][0]; y2=xyz_list[1][1]; - x3=xyz_list[2][0]; y3=xyz_list[2][1]; - - if(which_one==1){ - /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ - return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; - } - else if(which_one==2){ - /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ - return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; - } - else if(which_one==3){ - /*Get third area coordinate 1-area1-area2: */ - return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; - } - else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!")); -} -/*}}}*/ + + /* Recover input data: */ + if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); + if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"missing velocity_obs input parameter"); + } + if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"missing velocity input parameter"); + } + + /*Initialize velocities: */ + for(i=0;imeanvel/(obs_vx_list[i]+numpar->epsvel),(double)2); + scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),(double)2); + if(obs_vx_list[i]==0)scalex=0; + if(obs_vy_list[i]==0)scaley=0; + relative_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2)); + } + } + else if(fit==2){ + /*We are using a logarithmic misfit: */ + for (i=0;iepsvel; //epsvel to avoid velocity being nil. + obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. + logarithmic_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); + } + } + else{ + /*Not supported yet! : */ + throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit)); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + +#ifdef _ISSM_DEBUG_ + for (i=0;icontrol_type,"drag")==0){ + if (!shelf){ + + //noise dampening + GetParameterDerivativeValue(&dk[0], &k[0],&xyz_list[0][0], gauss_l1l2l3); + Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; + + } + } + else if (strcmp(numpar->control_type,"B")==0){ + if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ + throw ErrorException(__FUNCT__,"parameter B not found in input"); + } + //noise dampening + GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); + Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; + + //min dampening + GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3); + if(B_gausscm_mindmp_value){ + Jelem+=numpar->cm_mindmp_slope*B_gauss*Jdet*gauss_weight; + } + + //max dampening + if(B_gauss>numpar->cm_maxdmp_value){ + Jelem+=numpar->cm_maxdmp_slope*B_gauss*Jdet*gauss_weight; + } + } + else{ + throw ErrorException(__FUNCT__,exprintf("%s%s","unsupported control type: ",numpar->control_type)); + } + + /*Differents misfits are allowed: */ + if(fit==0){ + /*Compute absolute misfit at gaussian point: */ + GetParameterValue(&absolute, &absolute_list[0],gauss_l1l2l3); + + /*compute Misfit*/ + Jelem+=absolute*Jdet*gauss_weight; + } + else if(fit==1){ + /*Compute relative misfit at gaussian point: */ + GetParameterValue(&relative, &relative_list[0],gauss_l1l2l3); + + /*compute Misfit*/ + Jelem+=relative*Jdet*gauss_weight; + } + else if(fit==2){ + /*Compute logarithmic misfit at gaussian point: */ + GetParameterValue(&logarithmic, &logarithmic_list[0],gauss_l1l2l3); + + /*compute Misfit*/ + Jelem+=logarithmic*Jdet*gauss_weight; + } + else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!")); + + } +cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + /*Return: */ + return Jelem; +} +/*}}}*/ +/*FUNCTION MyRank {{{1*/ +int Tria::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ +/*FUNCTION NodeConfiguration {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::NodeConfiguration" +void Tria::NodeConfiguration(int* tria_node_ids,Node* tria_nodes[3],int* tria_node_offsets){ + + int i; + for(i=0;i<3;i++){ + node_ids[i]=tria_node_ids[i]; + nodes[i]=tria_nodes[i]; + node_offsets[i]=tria_node_offsets[i]; + } + +} +/*}}}*/ +/*FUNCTION NumparConfiguration {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::NumparConfiguration" +void Tria::NumparConfiguration(Numpar* tria_numpar,int tria_numpar_offset){ + + numpar=tria_numpar; + numpar_offset=tria_numpar_offset; + +} +/*}}}*/ +/*FUNCTION SurfaceNormal{{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::SurfaceNormal" + +void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ + + int i; + double v13[3]; + double v23[3]; + double normal[3]; + double normal_norm; + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0][i]-xyz_list[2][i]; + v23[i]=xyz_list[1][i]-xyz_list[2][i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + + normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) ); + + *(surface_normal)=normal[0]/normal_norm; + *(surface_normal+1)=normal[1]/normal_norm; + *(surface_normal+2)=normal[2]/normal_norm; + +} +/*}}}*/ +/*FUNCTION UpdateFromInputs {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::UpdateFromInputs" +void Tria::UpdateFromInputs(void* vinputs){ + + int dofs[1]={0}; + double temperature_list[3]; + double temperature_average; + double B_list[3]; + double B_average; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*Update internal data if inputs holds new values: */ + inputs->Recover("thickness",&h[0],1,dofs,3,(void**)nodes); + inputs->Recover("surface",&s[0],1,dofs,3,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,3,(void**)nodes); + inputs->Recover("drag",&k[0],1,dofs,3,(void**)nodes); + inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes); + inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes); + inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes); + + //Update material if necessary + if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){ + temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2])/3.0; + B_average=Paterson(temperature_average); + matice->SetB(B_average); + } + + if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){ + B_average=(B_list[0]+B_list[1]+B_list[2])/3.0; + matice->SetB(B_average); + } + +} +/*}}}*/