Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6157) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6158) @@ -540,11 +540,11 @@ void Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){ + bool already=false; int i,j; + int partition[NUMVERTICES]; + int offset[NUMVERTICES]; double area; double mean; - int partition[NUMVERTICES]; - int offset[NUMVERTICES]; double values[3]; - bool already=false; /*First, get the area: */ @@ -900,16 +900,13 @@ /*Intermediaries*/ int i,ig; - double Jdet; - double viscosity_complement; - double vx,vy,lambda,mu,thickness; - double cm_noisedmp; int doflist[NUMVERTICES]; + double vx,vy,lambda,mu,thickness,Jdet; + double cm_noisedmp,viscosity_complement; double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; double xyz_list[NUMVERTICES][3]; - double basis[3]; + double basis[3],epsilon[3]; double dbasis[NDOF2][NUMVERTICES]; + double grad_g[NUMVERTICES]; double grad[NUMVERTICES]={0.0}; - double grad_g[NUMVERTICES]; - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ GaussTria *gauss = NULL; @@ -954,12 +951,8 @@ } /*Add regularization term*/ - for (i=0;iweight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); - } - + for (i=0;iweight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - /*retrieve some parameters: */ + /*retrieve some parameters ands return if iceshelf: */ this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - - /*Get out if shelf*/ if(IsOnShelf())return; - - /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); GetDofList1(&doflist1[0]); @@ -1024,29 +1005,17 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); + /*Build alpha_complement_list: */ if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); else alpha_complement=0; - /*Recover alpha_complement and k: */ dragcoefficient_input->GetParameterValue(&drag, gauss); - - /*recover lambda and mu: */ adjointx_input->GetParameterValue(&lambda, gauss); adjointy_input->GetParameterValue(&mu, gauss); - - /*recover vx and vy: */ vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss); - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); - - /*Get k derivative: dk/dx */ dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); @@ -1065,5 +1034,4 @@ } - /*Add grade_g to global vector gradient: */ VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); @@ -1071,5 +1039,4 @@ delete gauss; delete friction; - } /*}}}*/ @@ -1082,5 +1049,4 @@ double gradient_g[NUMVERTICES]; - /*Retrieve dof list*/ GetDofList1(&doflist1[0]); @@ -1089,5 +1055,4 @@ for(int i=0;iGetInput(AdjointEnum); ISSMASSERT(adjoint_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); - thickness_input->GetParameterValue(&thickness, gauss); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(l1l2l3, gauss); - - for(i=0;iweight*l1l2l3[i]; - } - - VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -/*FUNCTION Tria::GradjVy{{{1*/ -void Tria::GradjVy(Vec gradient){ /*Intermediaries*/ @@ -1145,6 +1068,6 @@ double Dlambda[2]; double xyz_list[NUMVERTICES][3]; + double grade_g[NUMVERTICES] = {0.0}; GaussTria *gauss = NULL; - double grade_g[NUMVERTICES] = {0.0}; /* Get node coordinates and dof list: */ @@ -1162,7 +1085,49 @@ gauss->GaussPoint(ig); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); thickness_input->GetParameterValue(&thickness, gauss); + for(i=0;iweight*l1l2l3[i]; + } + + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjVy{{{1*/ +void Tria::GradjVy(Vec gradient){ + + /*Intermediaries*/ + int i,ig; + int doflist1[NUMVERTICES]; + double thickness,Jdet; + double l1l2l3[3]; + double Dlambda[2]; + double xyz_list[NUMVERTICES][3]; + double grade_g[NUMVERTICES] = {0.0}; + GaussTria *gauss = NULL; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + + /*Retrieve all inputs we will be needing: */ + Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); + thickness_input->GetParameterValue(&thickness, gauss); + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctions(l1l2l3, gauss); @@ -1181,7 +1146,7 @@ /*Intermediary*/ + int control_type; Input* input=NULL; double cm_min,cm_max; - int control_type; /*retrieve some parameters: */ @@ -1249,5 +1214,4 @@ ISSMERROR("control type %s not implemented yet",EnumToString(control_type)); } - } /*}}}*/ @@ -1255,8 +1219,8 @@ bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ - int i; + bool converged=true; + int i; 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 @@ -1280,5 +1244,4 @@ xfree((void**)&old_inputs); return converged; - } /*}}}*/ @@ -1310,5 +1273,4 @@ else ISSMERROR("object %s not supported yet",EnumToString(object_enum)); - } /*}}}*/ @@ -1354,10 +1316,7 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==RheologyBbarEnum){ - input=this->matice->inputs->GetInput(enum_type); - } - else{ - input=this->inputs->GetInput(enum_type); - } + if (enum_type==RheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); + else input=this->inputs->GetInput(enum_type); + if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumToString(enum_type)); @@ -1365,5 +1324,4 @@ * object out of the input, with the additional step and time information: */ this->results->AddObject((Object*)input->SpawnResult(step,time)); - } /*}}}*/ @@ -1371,17 +1329,15 @@ double Tria::MassFlux( double* segment,bool process_units){ - int i; - const int numdofs=2; - double mass_flux=0; - double xyz_list[NUMVERTICES][3]; + + int i; + double mass_flux=0; + double xyz_list[NUMVERTICES][3]; + double normal[2]; + double length,rho_ice; + double x1,y1,x2,y2,h1,h2; + double vx1,vx2,vy1,vy2; GaussTria* gauss_1=NULL; GaussTria* gauss_2=NULL; - double normal[2]; - double length; - double x1,y1,x2,y2; - double h1,h2; - double vx1,vx2,vy1,vy2; - double rho_ice; /*Get material parameters :*/ @@ -1403,21 +1359,17 @@ gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); - /*get normal of segment: */ normal[0]=cos(atan2(x1-x2,y2-y1)); normal[1]=sin(atan2(x1-x2,y2-y1)); - /*get length of segment: */ length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); - /*get thickness and velocity at two segment extremities: */ Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + thickness_input->GetParameterValue(&h1, gauss_1); thickness_input->GetParameterValue(&h2, gauss_2); - - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); vx_input->GetParameterValue(&vx1,gauss_1); vx_input->GetParameterValue(&vx2,gauss_2); - - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); vy_input->GetParameterValue(&vy1,gauss_1); vy_input->GetParameterValue(&vy2,gauss_2); @@ -1431,8 +1383,7 @@ mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum,this->parameters); - /*clean up*/ + /*clean up and return:*/ delete gauss_1; delete gauss_2; - return mass_flux; } @@ -1488,5 +1439,4 @@ /*Assign output pointers:*/ *pmaxvel=maxvel; - } /*}}}*/ @@ -1502,5 +1452,4 @@ /*Assign output pointers:*/ *pmaxvx=maxvx; - } /*}}}*/ @@ -1530,5 +1479,4 @@ /*Assign output pointers:*/ *pmaxvz=maxvz; - } /*}}}*/ @@ -1544,5 +1492,4 @@ /*Assign output pointers:*/ *pminvel=minvel; - } /*}}}*/ @@ -1558,5 +1505,4 @@ /*Assign output pointers:*/ *pminvx=minvx; - } /*}}}*/ @@ -1572,5 +1518,4 @@ /*Assign output pointers:*/ *pminvy=minvy; - } /*}}}*/ @@ -1586,5 +1531,4 @@ /*Assign output pointers:*/ *pminvz=minvz; - } /*}}}*/ @@ -1593,15 +1537,11 @@ /*intermediary: */ - double C; - double maxabsvx; - double maxabsvy; + int i; + double C,dt; + double dx,dy; double maxx,minx; double maxy,miny; + double maxabsvx,maxabsvy; double xyz_list[NUMVERTICES][3]; - double dx,dy; - int i; - - /*output: */ - double dt; /*get CFL coefficient:*/ @@ -1633,5 +1573,4 @@ return dt; - } /*}}}*/ @@ -1687,33 +1626,19 @@ double Tria::SurfaceAbsVelMisfit(bool process_units){ - int i; - - /* output: */ - double Jelem=0; - - /* node data: */ - const int numdof=2*NUMVERTICES; - double xyz_list[NUMVERTICES][3]; - - /* grid data: */ - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double misfit_square_list[NUMVERTICES]; - double misfit_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - - /* gaussian points: */ - int ig; + const int numdof=NDOF2*NUMVERTICES; + + int i,ig; + double Jelem=0; + double velocity_mag,obs_velocity_mag; + double meanvel, epsvel,misfit,Jdet; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double misfit_square_list[NUMVERTICES]; + double misfit_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; GaussTria *gauss=NULL; - - /* parameters: */ - double velocity_mag,obs_velocity_mag; - double misfit; - - /* Jacobian: */ - double Jdet; - double meanvel, epsvel; /*If on water, return 0: */ @@ -1758,7 +1683,5 @@ /*Apply weights to misfits*/ - for (i=0;iGetVertexId(); //vertices id start at column 3 of the patch. @@ -2220,6 +2063,4 @@ int i; - - /*output: */ int numrows = 0; int numnodes = 0; @@ -2238,5 +2079,4 @@ *pnumvertices=NUMVERTICES; *pnumnodes=numnodes; - } /*}}}*/ @@ -2250,5 +2090,4 @@ elementresult->ProcessUnits(this->parameters); } - } /*}}}*/ @@ -2256,19 +2095,13 @@ double Tria::SurfaceArea(void){ - int i; - - /* output: */ + int i; double S; - - /* node data: */ + double normal[3]; + double v13[3],v23[3]; double xyz_list[NUMVERTICES][3]; - double v13[3]; - double v23[3]; - double normal[3]; /*If on water, return 0: */ if(IsOnWater())return 0; - /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); @@ -2295,6 +2128,6 @@ int tria_node_ids[3]; int tria_vertex_ids[3]; + int tria_type; double nodeinputs[3]; - int tria_type; /*Checks if debuging*/ @@ -2504,10 +2337,9 @@ /*If sheet: surface = bed + thickness*/ else{ - + /*The bed does not change, update surface only s = b + h*/ this->inputs->DuplicateInput(BedEnum,SurfaceEnum); //1: copy bed into surface this->inputs->AXPY(SurfaceEnum,1.0,ThicknessEnum); //2: surface = surface + 1 * thickness } - } /*}}}*/ @@ -4632,16 +4464,10 @@ int i,j; + int count=0; int numberofdofs=0; - int count=0; - - /*output: */ int* doflist=NULL; - /*First, figure out size of doflist: */ - for(i=0;i<3;i++){ - numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - } - - /*Allocate: */ + /*First, figure out size of doflist and create it: */ + for(i=0;i<3;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); doflist=(int*)xmalloc(numberofdofs*sizeof(int)); @@ -4655,5 +4481,4 @@ /*Assign output pointers:*/ *pdoflist=doflist; - } /*}}}*/ @@ -4662,7 +4487,5 @@ int i; - for(i=0;i<3;i++){ - doflist[i]=nodes[i]->GetDofList1(); - } + for(i=0;i<3;i++) doflist[i]=nodes[i]->GetDofList1(); } @@ -4696,10 +4519,7 @@ void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){ - /*Intermediaries*/ double value[NUMVERTICES]; - GaussTria *gauss = NULL; - - /*Recover input*/ - Input* input=inputs->GetInput(enumtype); + GaussTria *gauss = NULL; + Input *input = inputs->GetInput(enumtype); /*Checks in debugging mode*/ @@ -4738,11 +4558,10 @@ void Tria::GetSolutionFromInputsDiagnosticHoriz(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]; - double vx; - double vy; GaussTria* gauss=NULL; @@ -4764,9 +4583,8 @@ vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - } - - /*Add value to global vector*/ + values[i*NDOF2+0]=vx; + values[i*NDOF2+1]=vy; + } + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); @@ -4774,5 +4592,4 @@ delete gauss; xfree((void**)&doflist); - } /*}}}*/ @@ -4780,12 +4597,10 @@ void Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){ - int i; - const int numdofpervertex=2; - const int numdof=numdofpervertex*NUMVERTICES; + const int numdof=NDOF2*NUMVERTICES; + + int i,dummy; int* doflist=NULL; + double vx,vy; double values[numdof]; - double vx; - double vy; - int dummy; GaussTria* gauss=NULL; @@ -4807,9 +4622,8 @@ vx_input->GetParameterValue(&vx,gauss); vy_input->GetParameterValue(&vy,gauss); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - } - - /*Add value to global vector*/ + values[i*NDOF2+0]=vx; + values[i*NDOF2+1]=vy; + } + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); @@ -4817,5 +4631,4 @@ delete gauss; xfree((void**)&doflist); - } /*}}}*/ @@ -4823,10 +4636,7 @@ void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){ /*Compute the 2d Strain Rate (3 components): - * - * epsilon=[exx eyy exy] - */ + * epsilon=[exx eyy exy] */ int i; - double epsilonvx[3]; double epsilonvy[3]; @@ -4843,5 +4653,4 @@ /*Sum all contributions*/ for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - } /*}}}*/ @@ -4849,34 +4658,20 @@ void Tria::GradjDragStokes(Vec gradient){ - int i,ig; - - /* node data: */ - double xyz_list[NUMVERTICES][3]; - int doflist1[NUMVERTICES]; - double dh1dh3[NDOF2][NUMVERTICES]; - - /* grid data: */ - double drag; - double alpha_complement; - Friction* friction=NULL; + int i,ig; + int drag_type,analysis_type; + int doflist1[NUMVERTICES]; + double bed,thickness,Neff; + double lambda,mu,xi,Jdet,vx,vy,vz; + double alpha_complement,drag,cm_noisedmp; + double surface_normal[3],bed_normal[3]; + double xyz_list[NUMVERTICES][3]; + double dh1dh3[NDOF2][NUMVERTICES]; + double dk[NDOF2]; + double l1l2l3[3]; + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + double grade_g[NUMVERTICES]={0.0}; + double grade_g_gaussian[NUMVERTICES]; + Friction* friction=NULL; GaussTria* gauss=NULL; - - /* parameters: */ - double vx,vy,vz; - double lambda,mu,xi; - double bed,thickness,Neff; - double surface_normal[3]; - double bed_normal[3]; - double dk[NDOF2]; - - /*element vector at the gaussian points: */ - double grade_g[NUMVERTICES]={0.0}; - double grade_g_gaussian[NUMVERTICES]; - double Jdet; - double l1l2l3[3]; - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - int drag_type; - double cm_noisedmp; - int analysis_type; /*retrive parameters: */ @@ -4964,10 +4759,8 @@ } - /*Add grade_g to global vector gradient: */ VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); delete friction; delete gauss; - } /*}}}*/ @@ -4975,11 +4768,10 @@ void Tria::InputUpdateFromSolutionAdjointBalancedthickness(double* solution){ - int i; - - const int numdofpervertex=1; - const int numdof=numdofpervertex*NUMVERTICES; - int* doflist=NULL; - double values[numdof]; - double lambda[NUMVERTICES]; + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double values[numdof]; + double lambda[NUMVERTICES]; /*Get dof list: */ @@ -4987,12 +4779,8 @@ /*Use the dof list to index into the solution vector: */ - for(i=0;iGetG(); GetParameterListOnVertices(&thickness[0],ThicknessEnum); - - for(i=0;i