Changeset 8646
- Timestamp:
- 06/16/11 14:13:46 (14 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r8643 r8646 3876 3876 /*FUNCTION Penta::Gradj {{{1*/ 3877 3877 void Penta::Gradj(Vec gradient,int control_type){ 3878 /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ 3878 3879 3879 3880 int i,approximation; … … 3883 3884 if(IsOnWater())return; 3884 3885 3886 /*First deal with ∂/∂alpha(KU-F)*/ 3885 3887 switch(control_type){ 3886 3888 … … 3962 3964 default: 3963 3965 _error_("control type %s not supported yet: ",EnumToStringx(control_type)); 3966 } 3967 3968 /*Now deal with ∂J/∂alpha*/ 3969 int *responses = NULL; 3970 int num_responses,resp; 3971 this->parameters->FindParam(&num_responses,NumResponsesEnum); 3972 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 3973 3974 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 3975 3976 case ThicknessAbsMisfitEnum: 3977 case SurfaceAbsVelMisfitEnum: 3978 case SurfaceRelVelMisfitEnum: 3979 case SurfaceLogVelMisfitEnum: 3980 case SurfaceLogVxVyMisfitEnum: 3981 case SurfaceAverageVelMisfitEnum: 3982 /*Nothing, J does not depends on the parameter being inverted for*/ 3983 break; 3984 case DragCoefficientAbsGradientEnum: 3985 if (!IsOnBed()) return; 3986 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 3987 tria->GradjDragGradient(gradient,resp); 3988 delete tria->matice; delete tria; 3989 break; 3990 default: 3991 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 3964 3992 } 3965 3993 } … … 5046 5074 /*Create inputs and add to DataSetInput*/ 5047 5075 DatasetInput* datasetinput=new DatasetInput(WeightsEnum); 5048 for( j=0;j<iomodel->num_cm_responses;j++){5049 for( i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1];5076 for(i=0;i<iomodel->num_cm_responses;i++){ 5077 for(j=0;j<6;j++)nodeinputs[j]=iomodel->weights[(penta_vertex_ids[j]-1)*iomodel->num_cm_responses+i]; 5050 5078 datasetinput->inputs->AddObject(new PentaVertexInput(WeightsEnum,nodeinputs)); 5051 5079 } -
issm/trunk/src/c/objects/Elements/Tria.cpp
r8643 r8646 1990 1990 } 1991 1991 break; 1992 case DragCoefficientAbsGradientEnum: 1993 /*Nothing in P vector*/ 1994 break; 1995 case ThicknessAbsGradientEnum: 1996 /*Nothing in P vector*/ 1997 break; 1998 case RheologyBbarAbsGradientEnum: 1999 /*Nothing in P vector*/ 2000 break; 1992 2001 default: 1993 2002 _error_("response %s not supported yet",EnumToStringx(responses[resp])); … … 2876 2885 /*FUNCTION Tria::Gradj {{{1*/ 2877 2886 void Tria::Gradj(Vec gradient,int control_type){ 2887 /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/ 2878 2888 2879 2889 /*If on water, grad = 0: */ 2880 2890 if(IsOnWater()) return; 2881 2891 2892 /*First deal with ∂/∂alpha(KU-F)*/ 2882 2893 switch(control_type){ 2883 2894 case DragCoefficientEnum: … … 2899 2910 _error_("%s%i","control type not supported yet: ",control_type); 2900 2911 } 2912 2913 /*Now deal with ∂J/∂alpha*/ 2914 int *responses = NULL; 2915 int num_responses,resp; 2916 this->parameters->FindParam(&num_responses,NumResponsesEnum); 2917 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 2918 2919 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 2920 2921 case ThicknessAbsMisfitEnum: 2922 case SurfaceAbsVelMisfitEnum: 2923 case SurfaceRelVelMisfitEnum: 2924 case SurfaceLogVelMisfitEnum: 2925 case SurfaceLogVxVyMisfitEnum: 2926 case SurfaceAverageVelMisfitEnum: 2927 /*Nothing, J does not depends on the parameter being inverted for*/ 2928 break; 2929 case DragCoefficientAbsGradientEnum: 2930 GradjDragGradient(gradient,resp); 2931 break; 2932 default: 2933 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 2934 } 2901 2935 } 2902 2936 /*}}}*/ … … 3157 3191 3158 3192 delete friction; 3193 delete gauss; 3194 } 3195 /*}}}*/ 3196 /*FUNCTION Tria::GradjDragGradient{{{1*/ 3197 void Tria::GradjDragGradient(Vec gradient, int weight_index){ 3198 3199 int i,ig; 3200 int doflist1[NUMVERTICES]; 3201 double Jdet,weight; 3202 double xyz_list[NUMVERTICES][3]; 3203 double dh1dh3[NDOF2][NUMVERTICES]; 3204 double dk[NDOF2]; 3205 double grade_g[NUMVERTICES]={0.0}; 3206 GaussTria *gauss=NULL; 3207 3208 /*Retrieve all inputs we will be needing: */ 3209 if(IsOnShelf())return; 3210 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 3211 GetDofList1(&doflist1[0]); 3212 Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); _assert_(dragcoefficient_input); 3213 Input* weights_input=inputs->GetInput(WeightsEnum); _assert_(weights_input); 3214 3215 /* Start looping on the number of gaussian points: */ 3216 gauss=new GaussTria(2); 3217 for (ig=gauss->begin();ig<gauss->end();ig++){ 3218 3219 gauss->GaussPoint(ig); 3220 3221 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 3222 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); 3223 weights_input->GetParameterValue(&weight,gauss,weight_index); 3224 3225 /*Build alpha_complement_list: */ 3226 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); 3227 3228 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ 3229 for (i=0;i<NUMVERTICES;i++){ 3230 grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); 3231 _assert_(!isnan(grade_g[i])); 3232 } 3233 } 3234 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 3235 3236 /*Clean up and return*/ 3159 3237 delete gauss; 3160 3238 } … … 3641 3719 /*Create inputs and add to DataSetInput*/ 3642 3720 DatasetInput* datasetinput=new DatasetInput(WeightsEnum); 3643 for( j=0;j<iomodel->num_cm_responses;j++){3644 for( i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1];3721 for(i=0;i<iomodel->num_cm_responses;i++){ 3722 for(j=0;j<3;j++)nodeinputs[j]=iomodel->weights[(tria_vertex_ids[j]-1)*iomodel->num_cm_responses+i]; 3645 3723 datasetinput->inputs->AddObject(new TriaVertexInput(WeightsEnum,nodeinputs)); 3646 3724 } -
issm/trunk/src/c/objects/Elements/Tria.h
r8643 r8646 93 93 void GradjDragMacAyeal(Vec gradient); 94 94 void GradjDragStokes(Vec gradient); 95 void GradjDragGradient(Vec gradient,int weight_index); 95 96 void GradjDhDtBalancedthickness(Vec gradient); 96 97 void GradjVxBalancedthickness(Vec gradient); -
issm/trunk/src/c/objects/Params/DoubleMatParam.cpp
r8224 r8646 152 152 153 153 /*DoubleMatParam virtual functions definitions: */ 154 /*FUNCTION DoubleMatParam::GetParameterValue {{{1*/154 /*FUNCTION DoubleMatParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){{{1*/ 155 155 void DoubleMatParam::GetParameterValue(double** pdoublearray,int* pM,int* pN){ 156 156 double* output=NULL; … … 163 163 if(pN) *pN=N; 164 164 *pdoublearray=output; 165 } 166 /*}}}*/ 167 /*FUNCTION DoubleMatParam::GetParameterValue(int** pintarray,int* pM,int* pN){{{1*/ 168 void DoubleMatParam::GetParameterValue(int** pintarray,int* pM,int* pN){ 169 #ifdef _SERIAL_ 170 int* output=NULL; 171 int i; 172 173 output=(int*)xmalloc((int)(M*N*sizeof(int))); 174 for(i=0;i<M*N;i++) output[i]=(int)value[i]; 175 176 /*Assign output pointers:*/ 177 if(pM) *pM=M; 178 if(pN) *pN=N; 179 *pintarray=output; 180 #else 181 _error_("DoubleMat of enum %i (%s) cannot return an array of int",enum_type,EnumToStringx(enum_type)); 182 #endif 165 183 } 166 184 /*}}}*/ -
issm/trunk/src/c/objects/Params/DoubleMatParam.h
r8600 r8646 54 54 void GetParameterValue(int* pinteger){_error_("DoubleMat param of enum %i (%s) cannot return an integer",enum_type,EnumToStringx(enum_type));} 55 55 void GetParameterValue(int** pintarray,int* pM){_error_("DoubleMat param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));} 56 void GetParameterValue(int** pintarray,int* pM,int* pN) {_error_("DoubleMat param of enum %i (%s) cannot return an array of integers",enum_type,EnumToStringx(enum_type));}56 void GetParameterValue(int** pintarray,int* pM,int* pN); 57 57 void GetParameterValue(double* pdouble){_error_("DoubleMat param of enum %i (%s) cannot return a double",enum_type,EnumToStringx(enum_type));} 58 58 void GetParameterValue(char** pstring){_error_("DoubleMat param of enum %i (%s) cannot return a string",enum_type,EnumToStringx(enum_type));} -
issm/trunk/src/c/objects/Params/DoubleVecParam.cpp
r8224 r8646 174 174 *pintarray=output; 175 175 #else 176 _error_("Double param of enum %i (%s) cannot return an array of double",enum_type,EnumToStringx(enum_type));176 _error_("DoubleVec param of enum %i (%s) cannot return an array of int",enum_type,EnumToStringx(enum_type)); 177 177 #endif 178 178 }
Note:
See TracChangeset
for help on using the changeset viewer.