Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6319) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6320) @@ -652,40 +652,22 @@ void Penta::ComputeBasalStress(Vec sigma_b){ - int i,j; - int doflist[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - double xyz_list_tria[3][3]; - - /*Parameters*/ - double rho_ice,gravity; - double bed_normal[3]; - double bed; - double basalforce[3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double stresstensor[6]={0.0}; - double viscosity; - int analysis_type; - - int dofv[3]={0,1,2}; - int dofp[1]={3}; - double Jdet2d; - - /*Gauss*/ - int ig; + int i,j,ig; + int dofv[3]={0,1,2}; + int dofp[1]={3}; + int analysis_type,approximation; + int doflist[NUMVERTICES]; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[3][3]; + double rho_ice,gravity,stokesreconditioning; + double pressure,viscosity,bed,Jdet2d; + double bed_normal[3]; + double basalforce[3]; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double stresstensor[6]={0.0}; + double sigma_xx,sigma_yy,sigma_zz; + double sigma_xy,sigma_xz,sigma_yz; + double surface=0,value=0; GaussPenta* gauss; - - /*Output*/ - double pressure; - double sigma_xx,sigma_yy,sigma_zz; - double sigma_xy,sigma_xz,sigma_yz; - double surface=0; - double value=0; - - /*flags: */ - int approximation; - - /*parameters: */ - double stokesreconditioning; /*retrive parameters: */ @@ -712,9 +694,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - for(i=0;i<3;i++){ - for(j=0;j<3;j++){ - xyz_list_tria[i][j]=xyz_list[i][j]; - } - } + for(i=0;i<3;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j]; /*Retrieve all inputs we will be needing: */ @@ -1114,7 +1092,7 @@ int i; + bool converged=true; Input** new_inputs=NULL; Input** old_inputs=NULL; - bool converged=true; new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs @@ -1140,5 +1118,4 @@ /*Return output*/ return converged; - } /*}}}*/ @@ -1189,7 +1166,8 @@ void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ - /*Intermediaries*/ int step,i; - + double xyz_list[NUMVERTICES][3]; + double Helem_list[NUMVERTICES]; + double zeros_list[NUMVERTICES]={0.0}; Penta* penta=NULL; Input* original_input=NULL; @@ -1199,8 +1177,4 @@ Input* total_thickness_input=NULL; Input* depth_averaged_input=NULL; - - double xyz_list[NUMVERTICES][3]; - double Helem_list[NUMVERTICES]; - double zeros_list[NUMVERTICES]={0.0}; /*recover parameters: */ @@ -1408,5 +1382,4 @@ /*Assign output pointers:*/ *pmaxvx=maxvx; - } /*}}}*/ @@ -1422,5 +1395,4 @@ /*Assign output pointers:*/ *pmaxvy=maxvy; - } /*}}}*/ @@ -1436,5 +1408,4 @@ /*Assign output pointers:*/ *pmaxvz=maxvz; - } /*}}}*/ @@ -1450,5 +1421,4 @@ /*Assign output pointers:*/ *pminvel=minvel; - } /*}}}*/ @@ -1464,5 +1434,4 @@ /*Assign output pointers:*/ *pminvx=minvx; - } /*}}}*/ @@ -1478,5 +1447,4 @@ /*Assign output pointers:*/ *pminvy=minvy; - } /*}}}*/ @@ -1492,5 +1460,4 @@ /*Assign output pointers:*/ *pminvz=minvz; - } /*}}}*/ @@ -1498,18 +1465,9 @@ double Penta::TimeAdapt(void){ - /*intermediary: */ - double C; - double maxabsvx; - double maxabsvy; - double maxabsvz; - double maxx,minx; - double maxy,miny; - double maxz,minz; - double dx,dy,dz; + int i; + double C,dx,dy,dz,dt; + double maxabsvx,maxabsvy,maxabsvz; + double maxx,minx,maxy,miny,maxz,minz; double xyz_list[NUMVERTICES][3]; - int i; - - /*output: */ - double dt; /*get CFL coefficient:*/ @@ -1547,14 +1505,11 @@ return dt; - } /*FUNCTION Penta::ThicknessAbsMisfit {{{1*/ double Penta::ThicknessAbsMisfit(bool process_units){ + int approximation; double J; - Tria* tria=NULL; - - /*inputs: */ - int approximation; + Tria* tria=NULL; /*retrieve inputs :*/ @@ -1574,9 +1529,7 @@ double Penta::SurfaceAbsVelMisfit(bool process_units){ + int approximation; double J; - Tria* tria=NULL; - - /*inputs: */ - int approximation; + Tria* tria=NULL; /*retrieve inputs :*/ @@ -1613,9 +1566,7 @@ double Penta::SurfaceRelVelMisfit(bool process_units){ + int approximation; double J; - Tria* tria=NULL; - - /*inputs: */ - int approximation; + Tria* tria=NULL; /*retrieve inputs :*/ @@ -1652,9 +1603,7 @@ double Penta::SurfaceLogVelMisfit(bool process_units){ + int approximation; double J; - Tria* tria=NULL; - - /*inputs: */ - int approximation; + Tria* tria=NULL; /*retrieve inputs :*/ @@ -1730,9 +1679,7 @@ double Penta::SurfaceAverageVelMisfit(bool process_units){ + int approximation; double J; - Tria* tria=NULL; - - /*inputs: */ - int approximation; + Tria* tria=NULL; /*retrieve inputs :*/ @@ -1769,8 +1716,6 @@ void Penta::PatchFill(int* pcount, Patch* patch){ - int i; - int count; + int i,count; int vertices_ids[6]; - /*recover pointer: */ @@ -1799,8 +1744,6 @@ int i; - - /*output: */ - int numrows = 0; - int numnodes = 0; + int numrows = 0; + int numnodes = 0; /*Go through all the results objects, and update the counters: */ @@ -1817,5 +1760,4 @@ *pnumvertices=NUMVERTICES; *pnumnodes=numnodes; - } /*}}}*/ @@ -1829,5 +1771,4 @@ elementresult->ProcessUnits(this->parameters); } - } /*}}}*/ @@ -1851,9 +1792,7 @@ double Penta::SurfaceArea(void){ + int approximation; double S; - Tria* tria=NULL; - - /*inputs: */ - int approximation; + Tria* tria=NULL; /*retrieve inputs :*/ @@ -1914,7 +1853,5 @@ /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<6;i++){ - penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab - } + for(i=0;i<6;i++) penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab /*Recover nodes ids needed to initialize the node hook.*/ @@ -2003,5 +1940,4 @@ break; } - } /*}}}*/ @@ -2010,5 +1946,5 @@ /*Intermediaries*/ - int i; + int i; double rho_ice,rho_water; @@ -2039,5 +1975,4 @@ this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum); //2: surface = surface + 1 * thickness } - } /*}}}*/ @@ -2048,6 +1983,5 @@ int i; - double v13[3]; - double v23[3]; + double v13[3],v23[3]; double normal[3]; double normal_norm; @@ -2061,5 +1995,4 @@ 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) ); @@ -2068,5 +2001,4 @@ *(bed_normal+1)=-normal[1]/normal_norm; *(bed_normal+2)=-normal[2]/normal_norm; - } /*}}}*/ @@ -3939,15 +3871,10 @@ void Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){ - int i,j; - int numberofdofs=0; - int count=0; - - /*output: */ + int i,j,count=0; + int numberofdofs=0; int* doflist=NULL; /*First, figure out size of doflist: */ - for(i=0;i<6;i++){ - numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - } + for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); /*Allocate: */ @@ -3963,5 +3890,4 @@ /*Assign output pointers:*/ *pdoflist=doflist; - } /*}}}*/ @@ -3970,8 +3896,5 @@ int i; - - for(i=0;i<6;i++){ - doflist[i]=nodes[i]->GetDofList1(); - } + for(i=0;i<6;i++) doflist[i]=nodes[i]->GetDofList1(); } @@ -3990,5 +3913,4 @@ /*return PentaRef field*/ return this->element_type; - } /*}}}*/ @@ -4101,14 +4023,12 @@ void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ - int i; - - const int numdofpervertex=2; - const int numdof=numdofpervertex*NUMVERTICES; + const int numdof=NDOF2*NUMVERTICES; + + int i; + int approximation; + int* doflist=NULL; + double vx,vy; + double values[numdof]; GaussPenta* gauss; - int* doflist=NULL; - double values[numdof]; - double vx; - double vy; - int approximation; /*Get approximation enum and dof list: */ @@ -4130,6 +4050,6 @@ vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; + values[i*NDOF2+0]=vx; + values[i*NDOF2+1]=vy; } @@ -4145,13 +4065,11 @@ void Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){ - int i; - - const int numdofpervertex=2; - const int numdof=numdofpervertex*NUMVERTICES; + const int numdof=NDOF2*NUMVERTICES; + + int i; + int* doflist=NULL; + double vx,vy; + double values[numdof]; GaussPenta* gauss=NULL; - int* doflist=NULL; - double values[numdof]; - double vx; - double vy; /*Get dof list: */ @@ -4164,11 +4082,10 @@ gauss=new GaussPenta(); for(i=0;iGaussVertex(i); vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; + values[i*NDOF2+0]=vx; + values[i*NDOF2+1]=vy; } @@ -4184,12 +4101,11 @@ void Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){ - int i; - - const int numdofpervertex=1; - const int numdof=numdofpervertex*NUMVERTICES; + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double vz; + double values[numdof]; GaussPenta* gauss=NULL; - int* doflist=NULL; - double values[numdof]; - double vz; /*Get dof list: */ @@ -4201,5 +4117,4 @@ gauss=new GaussPenta(); for(i=0;iGaussVertex(i); @@ -4219,13 +4134,12 @@ void Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){ - int i; - - const int numdofpervertex=4; - const int numdof=numdofpervertex*NUMVERTICES; - GaussPenta *gauss; + const int numdof=NDOF4*NUMVERTICES; + + int i; int* doflist=NULL; - double values[numdof]; double vx,vy,vz,p; double stokesreconditioning; + double values[numdof]; + GaussPenta *gauss; /*Get dof list: */ @@ -4248,8 +4162,8 @@ vz_input->GetParameterValue(&vz,gauss); p_input ->GetParameterValue(&p ,gauss); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - values[i*numdofpervertex+2]=vz; - values[i*numdofpervertex+3]=p/stokesreconditioning; + values[i*NDOF4+0]=vx; + values[i*NDOF4+1]=vy; + values[i*NDOF4+2]=vz; + values[i*NDOF4+3]=p/stokesreconditioning; } @@ -4265,13 +4179,11 @@ void Penta::GetSolutionFromInputsThermal(Vec solution){ - int i; - - const int numdofpervertex=1; - const int numdof=numdofpervertex*NUMVERTICES; - GaussPenta *gauss=NULL; + const int numdof=NDOF1*NUMVERTICES; + + int i; int* doflist=NULL; double values[numdof]; double temp; - + GaussPenta *gauss=NULL; /*Get dof list: */ @@ -4279,10 +4191,7 @@ Input* t_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(t_input); - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ gauss=new GaussPenta(); for(i=0;iGaussVertex(i); t_input->GetParameterValue(&temp,gauss); @@ -4327,5 +4236,4 @@ int i; - double epsilonvx[5]; double epsilonvy[5]; @@ -4342,5 +4250,4 @@ /*Sum all contributions*/ for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - } /*}}}*/ @@ -4353,5 +4260,4 @@ int i; - double epsilonvx[6]; double epsilonvy[6]; @@ -4370,5 +4276,4 @@ /*Sum all contributions*/ for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; - } /*}}}*/ @@ -4382,5 +4287,4 @@ penta=this; for(;;){ - /*Stop if we have reached the surface, else, take lower penta*/ if (penta->IsOnBed()) break; @@ -4399,7 +4303,8 @@ Penta* upper_penta=NULL; + upper_penta=(Penta*)neighbors[0]; //first one (0) under, second one (1) above + return upper_penta; - } /*}}}*/ @@ -4408,7 +4313,8 @@ Penta* upper_penta=NULL; + upper_penta=(Penta*)neighbors[1]; //first one under, second one above + return upper_penta; - } /*}}}*/ @@ -4417,12 +4323,10 @@ int i; + double z; double xyz_list[NUMVERTICES][3]; double z_list[NUMVERTICES]; - double z; GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - for(i=0;iGetParameterValue(&approximation,ApproximationEnum); /*If on water, skip: */ if(IsOnWater())return; + inputs->GetParameterValue(&approximation,ApproximationEnum); if (approximation==MacAyealApproximationEnum){ @@ -4476,5 +4378,4 @@ this->matice->inputs->DeleteInput(RheologyBbarEnum); } - } /*}}}*/ @@ -4482,24 +4383,16 @@ void Penta::GradjDrag(Vec gradient){ - int i; - Tria* tria=NULL; + int i,approximation; + double temp_gradient[6]={0,0,0,0,0,0}; + Tria* tria=NULL; TriaVertexInput* triavertexinput=NULL; - double temp_gradient[6]={0,0,0,0,0,0}; /*retrieve inputs :*/ - int approximation; inputs->GetParameterValue(&approximation,ApproximationEnum); - /*If on water, skip: */ - if(IsOnWater())return; - - /*If on shelf, skip: */ - if(IsOnShelf())return; - - /*Bail out if this element does not touch the bedrock: */ - if (!IsOnBed()) return; + /*If on water, on shelf or not on bed, skip: */ + if(IsOnWater()|| IsOnShelf() || !IsOnBed())return; if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum){ - /*MacAyeal or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. @@ -4508,5 +4401,4 @@ } else if (approximation==StokesApproximationEnum){ - /*Stokes*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. @@ -4518,6 +4410,4 @@ } else ISSMERROR("approximation %s not supported yet",EnumToString(approximation)); - - } /*}}}*/ @@ -4530,5 +4420,4 @@ /*Are we on the base, not on the surface?:*/ if(IsOnBed()){ - /*OK, we are on bed. we will follow the steps: * 1: find input and extrude it. @@ -4569,5 +4458,4 @@ /*Stop if we have reached the surface*/ if (penta->IsOnSurface()) break; - } } @@ -4579,5 +4467,4 @@ void Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ - /*Inputs*/ int approximation; @@ -4612,29 +4499,23 @@ } } - /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{1*/ void Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){ - int i; - - const int numdofpervertex=2; - const int numdof=numdofpervertex*NUMVERTICES; - int* doflist=NULL; - double values[numdof]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - int dummy; - double pressure[NUMVERTICES]; - double surface[NUMVERTICES]; - double rho_ice,g; - double xyz_list[NUMVERTICES][3]; - - double *vz_ptr = NULL; - Penta *penta = NULL; - - /*OK, we are on bed. Now recover results*/ + const int numdof=NDOF2*NUMVERTICES; + + int i,dummy; + double rho_ice,g; + double values[numdof]; + double vx[NUMVERTICES]; + double vy[NUMVERTICES]; + double vz[NUMVERTICES]; + double vel[NUMVERTICES]; + double pressure[NUMVERTICES]; + double surface[NUMVERTICES]; + double xyz_list[NUMVERTICES][3]; + int *doflist = NULL; + double *vz_ptr = NULL; + Penta *penta = NULL; /*Get dof list: */ @@ -4642,12 +4523,10 @@ /*Use the dof list to index into the solution vector: */ - for(i=0;iparameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - for(i=0;iinputs->GetParameterValue(&converged,ConvergedEnum); @@ -5441,8 +5269,8 @@ void Penta::InputUpdateFromSolutionOneDof(double* solution,int enum_type){ - const int numdofpervertex = 1; - const int numdof = numdofpervertex *NUMVERTICES; - int* doflist=NULL; - double values[numdof]; + const int numdof = NDOF1*NUMVERTICES; + + double values[numdof]; + int* doflist=NULL; /*Get dof list: */ @@ -5450,7 +5278,5 @@ /*Use the dof list to index into the solution vector: */ - for(int i=0;iinputs->AddInput(new PentaVertexInput(enum_type,values)); @@ -5538,5 +5363,4 @@ inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum); return onshelf; - } /*}}}*/ @@ -5547,5 +5371,4 @@ inputs->GetParameterValue(&onwater,ElementOnWaterEnum); return onwater; - } /*}}}*/ @@ -5556,5 +5379,4 @@ inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); return onsurface; - } /*}}}*/ @@ -5565,5 +5387,4 @@ inputs->GetParameterValue(&onbed,ElementOnBedEnum); return onbed; - } /*}}}*/ @@ -5574,6 +5395,6 @@ int i,node0,node1; int edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges + double length; double minlength=-1; - double length; for(i=0;i<9;i++){ @@ -5593,6 +5414,5 @@ void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){ - int i,j; - + int i,j; double Kii[24][24]; double Kib[24][3]; @@ -5630,9 +5450,5 @@ /*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]; - } - } + for(i=0;i<24;i++) for(j=0;j<24;j++) *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j]; } /*}}}*/ @@ -5640,6 +5456,5 @@ void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){ - int i,j; - + int i,j; double Pi[24]; double Pb[3]; @@ -5650,10 +5465,6 @@ /*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(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++){ @@ -5674,7 +5485,5 @@ /*Affect value to the reduced matrix */ - for(i=0;i<24;i++){ - *(Pe_reduced+i)=Pi[i]-Pright[i]; - } + for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i]; } /*}}}*/ @@ -5688,17 +5497,14 @@ Tria* Penta::SpawnTria(int g0, int g1, int g2){ - int i; - int analysis_counter; + int i,analysis_counter; + int indices[3]; + int zero=0; + Tria* tria = NULL; + Inputs* tria_inputs = NULL; + Results* tria_results = NULL; + Parameters* tria_parameters = NULL; /*go into parameters and get the analysis_counter: */ this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*out of grids g0,g1 and g2 from Penta, build a tria element: */ - Tria* tria=NULL; - int indices[3]; - int zero=0; - Parameters* tria_parameters=NULL; - Inputs* tria_inputs=NULL; - Results* tria_results=NULL; indices[0]=g0; @@ -5734,7 +5540,6 @@ void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ - int i; - double v13[3]; - double v23[3]; + int i; + double v13[3],v23[3]; double normal[3]; double normal_norm;