Index: /issm/trunk/src/c/objects/Elements/PentaRef.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5628) +++ /issm/trunk/src/c/objects/Elements/PentaRef.h (revision 5629) @@ -46,5 +46,7 @@ void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, double* gauss); void GetParameterValue(double* pvalue,double* plist,double* gauss); + void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, double* gauss); + void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");}; }; Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5628) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5629) @@ -29,4 +29,5 @@ this->parameters=NULL; this->results=NULL; + } /*}}}*/ @@ -34,6 +35,5 @@ Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels) :TriaRef(nummodels) - ,TriaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id - { //i is the element index + ,TriaHook(nummodels,index+1,iomodel->numberofelements+1){ /*id: */ this->id=tria_id; @@ -226,5 +226,9 @@ /*}}}*/ /*FUNCTION Tria::Id {{{1*/ -int Tria::Id(){ return id; } +int Tria::Id(){ + + return id; + +} /*}}}*/ /*FUNCTION Tria::Marshall {{{1*/ @@ -356,45 +360,45 @@ void Tria::InputUpdateFromSolution(double* solution){ + /*retrive parameters: */ int analysis_type; - - /*retrive parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHorizAnalysisEnum){ - InputUpdateFromSolutionDiagnosticHoriz( solution); - } - else if (analysis_type==DiagnosticHutterAnalysisEnum){ - InputUpdateFromSolutionDiagnosticHoriz( solution); - } - else if (analysis_type==AdjointHorizAnalysisEnum){ - InputUpdateFromSolutionAdjointHoriz( solution); - } - else if (analysis_type==BedSlopeXAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum); - } - else if (analysis_type==BedSlopeYAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum); - } - else if (analysis_type==SurfaceSlopeXAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum); - } - else if (analysis_type==SurfaceSlopeYAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum); - } - else if (analysis_type==PrognosticAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,ThicknessEnum); - } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,ThicknessEnum); - } - else if (analysis_type==AdjointBalancedthicknessAnalysisEnum){ - InputUpdateFromSolutionAdjointBalancedthickness( solution); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - InputUpdateFromSolutionOneDof(solution,VelEnum); - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + InputUpdateFromSolutionDiagnosticHoriz( solution); + break; + case DiagnosticHutterAnalysisEnum: + InputUpdateFromSolutionDiagnosticHoriz( solution); + break; + case AdjointHorizAnalysisEnum: + InputUpdateFromSolutionAdjointHoriz( solution); + break; + case BedSlopeXAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum); + break; + case BedSlopeYAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum); + break; + case SurfaceSlopeXAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum); + break; + case SurfaceSlopeYAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum); + break; + case PrognosticAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,ThicknessEnum); + break; + case BalancedthicknessAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,ThicknessEnum); + break; + case AdjointBalancedthicknessAnalysisEnum: + InputUpdateFromSolutionAdjointBalancedthickness( solution); + break; + case BalancedvelocitiesAnalysisEnum: + InputUpdateFromSolutionOneDof(solution,VelEnum); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); } } @@ -428,5 +432,4 @@ default: - ISSMERROR("type %i (%s) not implemented yet",type,EnumToString(type)); } @@ -446,6 +449,6 @@ void Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){ - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; int i,j; @@ -572,6 +575,6 @@ int i; - const int numgrids=3; - int doflist[numgrids]; + const int numvertices=3; + int doflist[numvertices]; double value; @@ -615,8 +618,8 @@ void Tria::ComputeStrainRate(Vec eps){ - int i; - const int numgrids=3; - int doflist[numgrids]; - double value; + int i; + const int numvertices = 3; + int doflist[numvertices]; + double value; /*plug local pressure values into global pressure vector: */ @@ -628,8 +631,7 @@ /*FUNCTION Tria::SetCurrentConfiguration {{{1*/ void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ - - int analysis_counter; /*go into parameters and get the analysis_counter: */ + int analysis_counter; parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); @@ -645,8 +647,7 @@ /*FUNCTION Tria::Configure {{{1*/ void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ - - int analysis_counter; /*go into parameters and get the analysis_counter: */ + int analysis_counter; parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); @@ -674,98 +675,66 @@ double Tria::RegularizeInversion(){ - int i; - - /* output: */ - double Jelem=0; - - /* node data: */ - const int numgrids=3; - const int numdof=2*numgrids; + /*constants: */ + const int numvertices=3; const int NDOF2=2; - double xyz_list[numgrids][3]; - - /* grid data: */ - 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 B_gauss; - double dhdt_gauss; - - /* parameters: */ - double dk[NDOF2]; - double dB[NDOF2]; - int control_type; - double cm_noisedmp; - - /* Jacobian: */ - double Jdet; - - /*inputs: */ - bool onwater; - bool shelf; - - /*retrieve inputs :*/ + const int numdof=NDOF2*numvertices; + + /* Intermediaries */ + int i,j,ig; + int control_type; + bool onwater,shelf; + double Jelem = 0; + double cm_noisedmp; + double Jdet; + double xyz_list[numvertices][3]; + double dk[NDOF2],dB[NDOF2]; + GaussTria *gauss = NULL; + + /*retrieve parameters and inputs*/ + this->parameters->FindParam(&control_type,ControlTypeEnum); + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); inputs->GetParameterValue(&onwater,ElementOnWaterEnum); inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + Input* drag_input=inputs->GetInput(DragCoefficientEnum); + Input* dhdt_input=inputs->GetInput(DhDtEnum); + Input* Bbar_input=matice->inputs->GetInput(RheologyBbarEnum); /*If on water, return 0: */ - if(onwater)return 0; - - /*retrieve some parameters: */ - this->parameters->FindParam(&control_type,ControlTypeEnum); - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussLegendreTria( &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; igbegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); /*Add Tikhonov regularization term to misfit*/ if (control_type==DragCoefficientEnum){ - if (shelf){ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); - Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; - - } + if (!shelf) return 0; //only shelf?... TO BE CHECKED + ISSMASSERT(drag_input); + drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; } else if (control_type==RheologyBbarEnum){ - matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBbarEnum); - Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; + ISSMASSERT(Bbar_input); + Bbar_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss); + Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; } else if (control_type==DhDtEnum){ - this->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],DhDtEnum); - Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; + ISSMASSERT(dhdt_input); + dhdt_input->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], gauss); + Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss->weight; } else{ ISSMERROR("%s%i","unsupported control type: ",control_type); } - - } - - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - - /*Return: */ + } + + /*clean-up and return: */ + delete gauss; return Jelem; } @@ -774,7 +743,6 @@ void Tria::CreateKMatrix(Mat Kgg){ + /*retrive parameters: */ int analysis_type; - - /*retrive parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -884,6 +852,6 @@ void Tria::GetBedList(double* bedlist){ - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; inputs->GetParameterValues(bedlist,&gaussgrids[0][0],3,BedEnum); @@ -953,6 +921,6 @@ void Tria::GetThicknessList(double* thickness_list){ - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],3,ThicknessEnum); @@ -1008,14 +976,14 @@ /* node data: */ - const int numgrids=3; + const int numvertices=3; const int NDOF1=1; const int NDOF2=2; - const int numdof=NDOF2*numgrids; - double xyz_list[numgrids][3]; - int doflist1[numgrids]; - double dh1dh3[NDOF2][numgrids]; + const int numdof=NDOF2*numvertices; + double xyz_list[numvertices][3]; + int doflist1[numvertices]; + double dh1dh3[NDOF2][numvertices]; /* grid data: */ - double B[numgrids]; + double B[numvertices]; @@ -1030,6 +998,6 @@ /*element vector at the gaussian points: */ - double grade_g[numgrids]={0.0}; - double grade_g_gaussian[numgrids]; + double grade_g[numvertices]={0.0}; + double grade_g_gaussian[numvertices]; /* Jacobian: */ @@ -1077,5 +1045,5 @@ /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); GetDofList1(&doflist1[0]); @@ -1128,5 +1096,5 @@ /*Build gradje_g_gaussian vector (actually -dJ/dB): */ - for (i=0;iGetAlphaComplement(&alpha2complement_list[i],&gauss[i][0],VxEnum,VyEnum); } @@ -1300,5 +1268,5 @@ /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iGetParameterValue(&lambda,&gauss[i][0]); gradient_g[i]=-lambda; @@ -1369,5 +1337,5 @@ /*Add grade_g to global vector gradient: */ - VecSetValues(gradient,numgrids,doflist1,(const double*)gradient_g,INSERT_VALUES); + VecSetValues(gradient,numvertices,doflist1,(const double*)gradient_g,INSERT_VALUES); } /*}}}*/ @@ -1568,8 +1536,8 @@ int i; - const int numgrids=3; + const int numvertices=3; const int numdofs=2; double mass_flux=0; - double xyz_list[numgrids][3]; + double xyz_list[numvertices][3]; double gauss_1[3]; double gauss_2[3]; @@ -1591,5 +1559,5 @@ /*Get xyz list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); /*get area coordinates of 0 and 1 locations: */ @@ -1630,7 +1598,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vx_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vx_values[numvertices]; double maxabsvx; @@ -1639,12 +1607,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numvertices,VxEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vx_values[0],numgrids,IuToExtEnum,VxEnum,this->parameters); + if(process_units)UnitConversion(&vx_values[0],numvertices,IuToExtEnum,VxEnum,this->parameters); /*now, compute maximum:*/ maxabsvx=fabs(vx_values[0]); - for(i=1;imaxabsvx)maxabsvx=fabs(vx_values[i]); } @@ -1659,7 +1627,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vy_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vy_values[numvertices]; double maxabsvy; @@ -1668,12 +1636,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numvertices,VyEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vy_values[0],numgrids,IuToExtEnum,VyEnum,this->parameters); + if(process_units)UnitConversion(&vy_values[0],numvertices,IuToExtEnum,VyEnum,this->parameters); /*now, compute maximum:*/ maxabsvy=fabs(vy_values[0]); - for(i=1;imaxabsvy)maxabsvy=fabs(vy_values[i]); } @@ -1688,7 +1656,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vz_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vz_values[numvertices]; double maxabsvz; @@ -1697,12 +1665,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numvertices,VzEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vz_values[0],numgrids,IuToExtEnum,VzEnum,this->parameters); + if(process_units)UnitConversion(&vz_values[0],numvertices,IuToExtEnum,VzEnum,this->parameters); /*now, compute maximum:*/ maxabsvz=fabs(vz_values[0]); - for(i=1;imaxabsvz)maxabsvz=fabs(vz_values[i]); } @@ -1717,7 +1685,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vel_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vel_values[numvertices]; double maxvel; @@ -1726,12 +1694,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numgrids,VelEnum); + inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numvertices,VelEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vel_values[0],numgrids,IuToExtEnum,VelEnum,this->parameters); + if(process_units)UnitConversion(&vel_values[0],numvertices,IuToExtEnum,VelEnum,this->parameters); /*now, compute maximum:*/ maxvel=vel_values[0]; - for(i=1;imaxvel)maxvel=vel_values[i]; } @@ -1747,7 +1715,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vx_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vx_values[numvertices]; double maxvx; @@ -1756,12 +1724,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numvertices,VxEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vx_values[0],numgrids,IuToExtEnum,VxEnum,this->parameters); + if(process_units)UnitConversion(&vx_values[0],numvertices,IuToExtEnum,VxEnum,this->parameters); /*now, compute maximum:*/ maxvx=vx_values[0]; - for(i=1;imaxvx)maxvx=vx_values[i]; } @@ -1777,7 +1745,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vy_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vy_values[numvertices]; double maxvy; @@ -1786,12 +1754,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numvertices,VyEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vy_values[0],numgrids,IuToExtEnum,VyEnum,this->parameters); + if(process_units)UnitConversion(&vy_values[0],numvertices,IuToExtEnum,VyEnum,this->parameters); /*now, compute maximum:*/ maxvy=vy_values[0]; - for(i=1;imaxvy)maxvy=vy_values[i]; } @@ -1807,7 +1775,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vz_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vz_values[numvertices]; double maxvz; @@ -1816,12 +1784,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numvertices,VzEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vz_values[0],numgrids,IuToExtEnum,VzEnum,this->parameters); + if(process_units)UnitConversion(&vz_values[0],numvertices,IuToExtEnum,VzEnum,this->parameters); /*now, compute maximum:*/ maxvz=vz_values[0]; - for(i=1;imaxvz)maxvz=vz_values[i]; } @@ -1837,7 +1805,7 @@ int i; int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; - double vel_values[numgrids]; + const int numvertices=3; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; + double vel_values[numvertices]; double minvel; @@ -1846,12 +1814,12 @@ /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numgrids,VelEnum); + inputs->GetParameterValues(&vel_values[0],&gaussgrids[0][0],numvertices,VelEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vel_values[0],numgrids,IuToExtEnum,VelEnum,this->parameters); + if(process_units)UnitConversion(&vel_values[0],numvertices,IuToExtEnum,VelEnum,this->parameters); /*now, compute minimum:*/ minvel=vel_values[0]; - for(i=1;iGetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numvertices,VxEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vx_values[0],numgrids,IuToExtEnum,VxEnum,this->parameters); + if(process_units)UnitConversion(&vx_values[0],numvertices,IuToExtEnum,VxEnum,this->parameters); /*now, compute minimum:*/ minvx=vx_values[0]; - for(i=1;iGetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numvertices,VyEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vy_values[0],numgrids,IuToExtEnum,VyEnum,this->parameters); + if(process_units)UnitConversion(&vy_values[0],numvertices,IuToExtEnum,VyEnum,this->parameters); /*now, compute minimum:*/ minvy=vy_values[0]; - for(i=1;iGetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numvertices,VzEnum); /*process units if requested: */ - if(process_units)UnitConversion(&vz_values[0],numgrids,IuToExtEnum,VzEnum,this->parameters); + if(process_units)UnitConversion(&vz_values[0],numvertices,IuToExtEnum,VzEnum,this->parameters); /*now, compute minimum:*/ minvz=vz_values[0]; - for(i=1;iparameters); + if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceAbsVelMisfitEnum,this->parameters); /*Apply weights to misfits*/ - for (i=0;iparameters); + if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceRelVelMisfitEnum,this->parameters); /*Apply weights to misfits*/ - for (i=0;iparameters); + if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceLogVelMisfitEnum,this->parameters); /*Apply weights to misfits*/ - for (i=0;iparameters); + if(process_units)UnitConversion(&misfit_list[0],numvertices,IuToExtEnum,SurfaceLogVxVyMisfitEnum,this->parameters); /*Apply weights to misfits*/ - for (i=0;iparameters); + if(process_units)UnitConversion(&misfit_square_list[0],numvertices,IuToExtEnum,SurfaceAverageVelMisfitEnum,this->parameters); /*Take the square root, and scale by surface: */ - for (i=0;iparameters->FindParam(&dim,DimEnum); @@ -3293,10 +3261,10 @@ /* node data: */ - const int numgrids=3; + const int numvertices=3; const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; + const int numdof=NDOF1*numvertices; + double xyz_list[numvertices][3]; int* doflist=NULL; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + double gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; /* gaussian points: */ @@ -3310,7 +3278,7 @@ /* matrices: */ - double L[numgrids]; - double B[2][numgrids]; - double Bprime[2][numgrids]; + double L[numvertices]; + double B[2][numvertices]; + double Bprime[2][numvertices]; double DL[2][2]={0.0}; double DLprime[2][2]={0.0}; @@ -3347,5 +3315,5 @@ /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices); GetDofList(&doflist); @@ -3363,5 +3331,5 @@ /*Modify z so that it reflects the surface*/ surface_input->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3); - for(i=0;iGetParameterValue(&weight, &gauss_l1l2l3[0]); - for (i=0;iGetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED } vx_input->GetParameterValues(&vx_list[0],&gauss[0][0],3); //TO BE DELETED vy_input->GetParameterValues(&vy_list[0],&gauss[0][0],3); //TO BE DELETED - for(i=0;i