Changeset 8607
- Timestamp:
- 06/10/11 16:54:02 (14 years ago)
- Location:
- issm/trunk/src
- Files:
-
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.cpp
r8224 r8607 11 11 #include "../Responsex/Responsex.h" 12 12 13 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters ,int response){13 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters){ 14 14 15 15 /*Intermediary*/ 16 16 int i; 17 double fit;17 int num_responses; 18 18 double S; 19 19 Element* element=NULL; 20 int* responses=NULL; 20 21 21 22 /*output: */ 22 double J; 23 double J,Jplus; 24 25 26 /*Recover parameters*/ 27 parameters->FindParam(&num_responses,NumResponsesEnum); 28 parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); 29 30 /*Get response*/ 31 J=0; 32 for(int i=0;i<num_responses;i++){ 33 Responsex(&Jplus,elements,nodes,vertices,loads,materials,parameters,EnumToStringx(responses[i]),false,i); //False means DO NOT process units 34 J+=Jplus; 35 } 36 37 /*REST TO BE DELETED*/ 38 39 /*Add Regularization terms: */ 23 40 double Jreg=0; 24 41 double Jreg_sum; 25 26 /*Get response*/27 Responsex(&J,elements,nodes,vertices,loads,materials,parameters,EnumToStringx(response),false); //False means DO NOT process units28 29 /*Add Regularization terms: */30 42 for (i=0;i<elements->Size();i++){ 31 43 element=(Element*)elements->GetObjectByOffset(i); … … 39 51 40 52 /*Assign output pointers: */ 53 xfree((void**)&responses); 41 54 *pJ=J; 42 55 } -
issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h
r5281 r8607 10 10 11 11 /* local prototypes: */ 12 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters ,int response);12 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters); 13 13 14 14 #endif -
issm/trunk/src/c/modules/DakotaResponsesx/DakotaResponsesx.cpp
r8303 r8607 76 76 77 77 //Responsex(responses_pointer,elements,nodes, vertices,loads,materials, parameters,root,process_units); 78 Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units );78 Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units,0);//0 is the index for weights 79 79 80 80 if(my_rank==0){ … … 95 95 96 96 /*perfectly normal response function: */ 97 Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units );97 Responsex(&femmodel_response,elements,nodes, vertices,loads,materials, parameters,root,process_units,0);//0 is the weight index 98 98 99 99 if(my_rank==0){ -
issm/trunk/src/c/modules/Responsex/Responsex.cpp
r8224 r8607 16 16 #include "../modules.h" 17 17 18 void Responsex(double* responses,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units ){18 void Responsex(double* responses,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units,int weight_index){ 19 19 20 20 switch (StringToEnumx(response_descriptor)){ 21 21 22 case MinVelEnum: MinVelx(responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;23 case MaxVelEnum: MaxVelx(responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;24 case MinVxEnum: MinVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);break;25 case MaxVxEnum: MaxVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);break;26 case MaxAbsVxEnum: MaxAbsVxx(responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;27 case MinVyEnum: MinVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);break;28 case MaxVyEnum: MaxVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);break;29 case MaxAbsVyEnum: MaxAbsVyx(responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;30 case MinVzEnum: MinVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);break;31 case MaxVzEnum: MaxVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units);break;32 case MaxAbsVzEnum: MaxAbsVzx(responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;33 case MassFluxEnum: MassFluxx(responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;34 case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;35 case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;36 case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;37 case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;38 case SurfaceAverageVelMisfitEnum: SurfaceAverageVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;39 case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break;22 case MinVelEnum: MinVelx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 23 case MaxVelEnum: MaxVelx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 24 case MinVxEnum: MinVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 25 case MaxVxEnum: MaxVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 26 case MaxAbsVxEnum: MaxAbsVxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 27 case MinVyEnum: MinVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 28 case MaxVyEnum: MaxVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 29 case MaxAbsVyEnum: MaxAbsVyx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 30 case MinVzEnum: MinVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 31 case MaxVzEnum: MaxVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 32 case MaxAbsVzEnum: MaxAbsVzx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 33 case MassFluxEnum: MassFluxx( responses, elements,nodes, vertices, loads, materials, parameters,process_units); break; 34 case SurfaceAbsVelMisfitEnum: SurfaceAbsVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 35 case SurfaceRelVelMisfitEnum: SurfaceRelVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 36 case SurfaceLogVelMisfitEnum: SurfaceLogVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 37 case SurfaceLogVxVyMisfitEnum: SurfaceLogVxVyMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 38 case SurfaceAverageVelMisfitEnum:SurfaceAverageVelMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 39 case ThicknessAbsMisfitEnum: ThicknessAbsMisfitx( responses, elements,nodes, vertices, loads, materials, parameters,process_units,weight_index); break; 40 40 default: _error_(" response descriptor \"%s\" not supported yet!",response_descriptor); break; 41 41 } -
issm/trunk/src/c/modules/Responsex/Responsex.h
r5473 r8607 9 9 #include "../../Container/Container.h" 10 10 11 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units );11 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,char* response_descriptor,bool process_units,int weight_index); 12 12 13 13 #endif /* _RESPONSESXX_H */ -
issm/trunk/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp
r5284 r8607 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 11 12 void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units ){12 void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){ 13 13 14 14 /*Intermediary*/ … … 23 23 for (i=0;i<elements->Size();i++){ 24 24 element=(Element*)elements->GetObjectByOffset(i); 25 J+=element->SurfaceAbsVelMisfit(process_units );25 J+=element->SurfaceAbsVelMisfit(process_units,weight_index); 26 26 } 27 27 -
issm/trunk/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h
r5284 r8607 10 10 11 11 /* local prototypes: */ 12 void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units );12 void SurfaceAbsVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); 13 13 14 14 #endif -
issm/trunk/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp
r5414 r8607 11 11 #include "../SurfaceAreax/SurfaceAreax.h" 12 12 13 void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units ){13 void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){ 14 14 15 15 /*Intermediary*/ … … 27 27 for (i=0;i<elements->Size();i++){ 28 28 element=(Element*)elements->GetObjectByOffset(i); 29 J+=element->SurfaceAverageVelMisfit(process_units );29 J+=element->SurfaceAverageVelMisfit(process_units,weight_index); 30 30 } 31 31 -
issm/trunk/src/c/modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h
r5284 r8607 10 10 11 11 /* local prototypes: */ 12 void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units );12 void SurfaceAverageVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); 13 13 14 14 #endif -
issm/trunk/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp
r5284 r8607 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 11 12 void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units ){12 void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){ 13 13 14 14 /*Intermediary*/ … … 23 23 for (i=0;i<elements->Size();i++){ 24 24 element=(Element*)elements->GetObjectByOffset(i); 25 J+=element->SurfaceLogVelMisfit(process_units );25 J+=element->SurfaceLogVelMisfit(process_units,weight_index); 26 26 } 27 27 -
issm/trunk/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h
r5284 r8607 10 10 11 11 /* local prototypes: */ 12 void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units );12 void SurfaceLogVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); 13 13 14 14 #endif -
issm/trunk/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp
r5284 r8607 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 11 12 void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units ){12 void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){ 13 13 14 14 /*Intermediary*/ … … 23 23 for (i=0;i<elements->Size();i++){ 24 24 element=(Element*)elements->GetObjectByOffset(i); 25 J+=element->SurfaceLogVxVyMisfit(process_units );25 J+=element->SurfaceLogVxVyMisfit(process_units,weight_index); 26 26 } 27 27 -
issm/trunk/src/c/modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h
r5284 r8607 10 10 11 11 /* local prototypes: */ 12 void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units );12 void SurfaceLogVxVyMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); 13 13 14 14 #endif -
issm/trunk/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp
r5284 r8607 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 11 12 void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units ){12 void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units,int weight_index){ 13 13 14 14 /*Intermediary*/ … … 23 23 for (i=0;i<elements->Size();i++){ 24 24 element=(Element*)elements->GetObjectByOffset(i); 25 J+=element->SurfaceRelVelMisfit(process_units );25 J+=element->SurfaceRelVelMisfit(process_units,weight_index); 26 26 } 27 27 -
issm/trunk/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h
r5284 r8607 10 10 11 11 /* local prototypes: */ 12 void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units );12 void SurfaceRelVelMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); 13 13 14 14 #endif -
issm/trunk/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp
r5286 r8607 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 11 12 void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units ){12 void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,bool process_units, int weight_index){ 13 13 14 14 /*Intermediary*/ … … 23 23 for (i=0;i<elements->Size();i++){ 24 24 element=(Element*)elements->GetObjectByOffset(i); 25 J+=element->ThicknessAbsMisfit(process_units );25 J+=element->ThicknessAbsMisfit(process_units,weight_index); 26 26 } 27 27 -
issm/trunk/src/c/modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h
r5286 r8607 10 10 11 11 /* local prototypes: */ 12 void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units );12 void ThicknessAbsMisfitx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,bool process_units,int weight_index); 13 13 14 14 #endif -
issm/trunk/src/c/objects/Elements/Element.h
r8365 r8607 43 43 virtual void GradjDrag(Vec gradient)=0; 44 44 virtual void GradjB(Vec gradient)=0; 45 virtual double ThicknessAbsMisfit(bool process_units )=0;46 virtual double SurfaceAbsVelMisfit(bool process_units )=0;47 virtual double SurfaceRelVelMisfit(bool process_units )=0;48 virtual double SurfaceLogVelMisfit(bool process_units )=0;49 virtual double SurfaceLogVxVyMisfit(bool process_units )=0;50 virtual double SurfaceAverageVelMisfit(bool process_units )=0;45 virtual double ThicknessAbsMisfit(bool process_units ,int weight_index)=0; 46 virtual double SurfaceAbsVelMisfit(bool process_units ,int weight_index)=0; 47 virtual double SurfaceRelVelMisfit(bool process_units ,int weight_index)=0; 48 virtual double SurfaceLogVelMisfit(bool process_units ,int weight_index)=0; 49 virtual double SurfaceLogVxVyMisfit(bool process_units,int weight_index)=0; 50 virtual double SurfaceAverageVelMisfit(bool process_units,int weight_index)=0; 51 51 virtual double RegularizeInversion(void)=0; 52 52 virtual double SurfaceArea(void)=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r8601 r8607 6845 6845 /*}}}*/ 6846 6846 /*FUNCTION Penta::SurfaceAverageVelMisfit {{{1*/ 6847 double Penta::SurfaceAverageVelMisfit(bool process_units ){6847 double Penta::SurfaceAverageVelMisfit(bool process_units,int weight_index){ 6848 6848 6849 6849 int approximation; … … 6868 6868 * and compute SurfaceAverageVelMisfit*/ 6869 6869 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 6870 J=tria->SurfaceAverageVelMisfit(process_units );6870 J=tria->SurfaceAverageVelMisfit(process_units,weight_index); 6871 6871 delete tria->matice; delete tria; 6872 6872 return J; … … 6875 6875 6876 6876 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 6877 J=tria->SurfaceAverageVelMisfit(process_units );6877 J=tria->SurfaceAverageVelMisfit(process_units,weight_index); 6878 6878 delete tria->matice; delete tria; 6879 6879 return J; … … 6882 6882 /*}}}*/ 6883 6883 /*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/ 6884 double Penta::SurfaceAbsVelMisfit(bool process_units ){6884 double Penta::SurfaceAbsVelMisfit(bool process_units,int weight_index){ 6885 6885 6886 6886 int approximation; … … 6905 6905 * and compute SurfaceAbsVelMisfit*/ 6906 6906 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 6907 J=tria->SurfaceAbsVelMisfit(process_units );6907 J=tria->SurfaceAbsVelMisfit(process_units,weight_index); 6908 6908 delete tria->matice; delete tria; 6909 6909 return J; … … 6912 6912 6913 6913 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 6914 J=tria->SurfaceAbsVelMisfit(process_units );6914 J=tria->SurfaceAbsVelMisfit(process_units,weight_index); 6915 6915 delete tria->matice; delete tria; 6916 6916 return J; … … 6919 6919 /*}}}*/ 6920 6920 /*FUNCTION Penta::SurfaceLogVelMisfit {{{1*/ 6921 double Penta::SurfaceLogVelMisfit(bool process_units ){6921 double Penta::SurfaceLogVelMisfit(bool process_units,int weight_index){ 6922 6922 6923 6923 int approximation; … … 6942 6942 * and compute SurfaceLogVelMisfit*/ 6943 6943 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 6944 J=tria->SurfaceLogVelMisfit(process_units );6944 J=tria->SurfaceLogVelMisfit(process_units,weight_index); 6945 6945 delete tria->matice; delete tria; 6946 6946 return J; … … 6949 6949 6950 6950 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 6951 J=tria->SurfaceLogVelMisfit(process_units );6951 J=tria->SurfaceLogVelMisfit(process_units,weight_index); 6952 6952 delete tria->matice; delete tria; 6953 6953 return J; … … 6956 6956 /*}}}*/ 6957 6957 /*FUNCTION Penta::SurfaceLogVxVyMisfit {{{1*/ 6958 double Penta::SurfaceLogVxVyMisfit(bool process_units ){6958 double Penta::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ 6959 6959 6960 6960 double J; … … 6981 6981 * and compute SurfaceLogVxVyMisfit*/ 6982 6982 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 6983 J=tria->SurfaceLogVxVyMisfit(process_units );6983 J=tria->SurfaceLogVxVyMisfit(process_units,weight_index); 6984 6984 delete tria->matice; delete tria; 6985 6985 return J; … … 6988 6988 6989 6989 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 6990 J=tria->SurfaceLogVxVyMisfit(process_units );6990 J=tria->SurfaceLogVxVyMisfit(process_units,weight_index); 6991 6991 delete tria->matice; delete tria; 6992 6992 return J; … … 7019 7019 /*}}}*/ 7020 7020 /*FUNCTION Penta::SurfaceRelVelMisfit {{{1*/ 7021 double Penta::SurfaceRelVelMisfit(bool process_units ){7021 double Penta::SurfaceRelVelMisfit(bool process_units,int weight_index){ 7022 7022 7023 7023 int approximation; … … 7042 7042 * and compute SurfaceRelVelMisfit*/ 7043 7043 tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face). 7044 J=tria->SurfaceRelVelMisfit(process_units );7044 J=tria->SurfaceRelVelMisfit(process_units,weight_index); 7045 7045 delete tria->matice; delete tria; 7046 7046 return J; … … 7049 7049 7050 7050 tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face). 7051 J=tria->SurfaceRelVelMisfit(process_units );7051 J=tria->SurfaceRelVelMisfit(process_units,weight_index); 7052 7052 delete tria->matice; delete tria; 7053 7053 return J; … … 7100 7100 } 7101 7101 /*FUNCTION Penta::ThicknessAbsMisfit {{{1*/ 7102 double Penta::ThicknessAbsMisfit(bool process_units ){7102 double Penta::ThicknessAbsMisfit(bool process_units,int weight_index){ 7103 7103 7104 7104 int approximation; … … 7113 7113 _error_("Not implemented yet"); 7114 7114 7115 tria=(Tria*)SpawnTria( 3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).7116 J=tria->ThicknessAbsMisfit(process_units );7115 tria=(Tria*)SpawnTria(0,1,2); 7116 J=tria->ThicknessAbsMisfit(process_units,weight_index); 7117 7117 delete tria->matice; delete tria; 7118 7118 return J; -
issm/trunk/src/c/objects/Elements/Penta.h
r8496 r8607 118 118 void MinVy(double* pminvy, bool process_units); 119 119 void MinVz(double* pminvz, bool process_units); 120 double ThicknessAbsMisfit( bool process_units);121 double SurfaceAbsVelMisfit( bool process_units);122 double SurfaceRelVelMisfit( bool process_units);123 double SurfaceLogVelMisfit( bool process_units);124 double SurfaceLogVxVyMisfit( bool process_units);125 double SurfaceAverageVelMisfit(bool process_units );120 double ThicknessAbsMisfit( bool process_units,int weight_index); 121 double SurfaceAbsVelMisfit( bool process_units,int weight_index); 122 double SurfaceRelVelMisfit( bool process_units,int weight_index); 123 double SurfaceLogVelMisfit( bool process_units,int weight_index); 124 double SurfaceLogVxVyMisfit( bool process_units,int weight_index); 125 double SurfaceAverageVelMisfit(bool process_units,int weight_index); 126 126 void PatchFill(int* pcount, Patch* patch); 127 127 void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r8603 r8607 1706 1706 1707 1707 /*Get all parameters at gaussian point*/ 1708 weights_input->GetParameterValue(&weight,gauss,0);1709 1708 vx_input->GetParameterValue(&vx,gauss); 1710 1709 vy_input->GetParameterValue(&vy,gauss); … … 1714 1713 1715 1714 /*Loop over all requested responses*/ 1716 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 1717 1718 case SurfaceAbsVelMisfitEnum: 1719 /* 1720 * 1 [ 2 2 ] 1721 * J = --- | (u - u ) + (v - v ) | 1722 * 2 [ obs obs ] 1723 * 1724 * dJ 1725 * DU = - -- = (u - u ) 1726 * du obs 1727 */ 1728 for (i=0;i<NUMVERTICES;i++){ 1729 dux=vxobs-vx; 1730 duy=vyobs-vy; 1731 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1732 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1733 } 1734 break; 1735 case SurfaceRelVelMisfitEnum: 1736 /* 1737 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1738 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1739 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1740 * obs obs 1741 * 1742 * dJ \bar{v}^2 1743 * DU = - -- = ------------- (u - u ) 1744 * du (u + eps)^2 obs 1745 * obs 1746 */ 1747 for (i=0;i<NUMVERTICES;i++){ 1748 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1749 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1750 dux=scalex*(vxobs-vx); 1751 duy=scaley*(vyobs-vy); 1752 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1753 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1754 } 1755 break; 1756 case SurfaceLogVelMisfitEnum: 1757 /* 1758 * [ vel + eps ] 2 1759 * J = 4 \bar{v}^2 | log ( ----------- ) | 1760 * [ vel + eps ] 1761 * obs 1762 * 1763 * dJ 2 * log(...) 1764 * DU = - -- = - 4 \bar{v}^2 ------------- u 1765 * du vel^2 + eps 1766 * 1767 */ 1768 for (i=0;i<NUMVERTICES;i++){ 1769 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1770 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1771 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1772 dux=scale*vx; 1773 duy=scale*vy; 1774 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1775 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1776 } 1777 break; 1778 case SurfaceAverageVelMisfitEnum: 1779 /* 1780 * 1 2 2 1781 * J = --- sqrt( (u - u ) + (v - v ) ) 1782 * S obs obs 1783 * 1784 * dJ 1 1 1785 * DU = - -- = - --- ----------- * 2 (u - u ) 1786 * du S 2 sqrt(...) obs 1787 */ 1788 for (i=0;i<NUMVERTICES;i++){ 1789 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1790 dux=scale*(vxobs-vx); 1791 duy=scale*(vyobs-vy); 1792 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1793 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1794 } 1795 break; 1796 case SurfaceLogVxVyMisfitEnum: 1797 /* 1798 * 1 [ |u| + eps 2 |v| + eps 2 ] 1799 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1800 * 2 [ |u |+ eps |v |+ eps ] 1801 * obs obs 1802 * dJ 1 u 1 1803 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1804 * du |u| + eps |u| u + eps 1805 */ 1806 for (i=0;i<NUMVERTICES;i++){ 1807 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1808 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1809 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1810 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1811 } 1812 break; 1813 default: 1814 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1715 for(resp=0;resp<num_responses;resp++){ 1716 1717 weights_input->GetParameterValue(&weight,gauss,resp); 1718 1719 switch(responses[resp]){ 1720 case SurfaceAbsVelMisfitEnum: 1721 /* 1722 * 1 [ 2 2 ] 1723 * J = --- | (u - u ) + (v - v ) | 1724 * 2 [ obs obs ] 1725 * 1726 * dJ 1727 * DU = - -- = (u - u ) 1728 * du obs 1729 */ 1730 for (i=0;i<NUMVERTICES;i++){ 1731 dux=vxobs-vx; 1732 duy=vyobs-vy; 1733 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1734 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1735 } 1736 break; 1737 case SurfaceRelVelMisfitEnum: 1738 /* 1739 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1740 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1741 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1742 * obs obs 1743 * 1744 * dJ \bar{v}^2 1745 * DU = - -- = ------------- (u - u ) 1746 * du (u + eps)^2 obs 1747 * obs 1748 */ 1749 for (i=0;i<NUMVERTICES;i++){ 1750 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1751 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1752 dux=scalex*(vxobs-vx); 1753 duy=scaley*(vyobs-vy); 1754 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1755 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1756 } 1757 break; 1758 case SurfaceLogVelMisfitEnum: 1759 /* 1760 * [ vel + eps ] 2 1761 * J = 4 \bar{v}^2 | log ( ----------- ) | 1762 * [ vel + eps ] 1763 * obs 1764 * 1765 * dJ 2 * log(...) 1766 * DU = - -- = - 4 \bar{v}^2 ------------- u 1767 * du vel^2 + eps 1768 * 1769 */ 1770 for (i=0;i<NUMVERTICES;i++){ 1771 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1772 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1773 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1774 dux=scale*vx; 1775 duy=scale*vy; 1776 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1777 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1778 } 1779 break; 1780 case SurfaceAverageVelMisfitEnum: 1781 /* 1782 * 1 2 2 1783 * J = --- sqrt( (u - u ) + (v - v ) ) 1784 * S obs obs 1785 * 1786 * dJ 1 1 1787 * DU = - -- = - --- ----------- * 2 (u - u ) 1788 * du S 2 sqrt(...) obs 1789 */ 1790 for (i=0;i<NUMVERTICES;i++){ 1791 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1792 dux=scale*(vxobs-vx); 1793 duy=scale*(vyobs-vy); 1794 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1795 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1796 } 1797 break; 1798 case SurfaceLogVxVyMisfitEnum: 1799 /* 1800 * 1 [ |u| + eps 2 |v| + eps 2 ] 1801 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1802 * 2 [ |u |+ eps |v |+ eps ] 1803 * obs obs 1804 * dJ 1 u 1 1805 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1806 * du |u| + eps |u| u + eps 1807 */ 1808 for (i=0;i<NUMVERTICES;i++){ 1809 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1810 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1811 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1812 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1813 } 1814 break; 1815 default: 1816 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1817 } 1815 1818 } 1816 1819 } … … 1870 1873 1871 1874 /*Get all parameters at gaussian point*/ 1872 weights_input->GetParameterValue(&weight,gauss,0);1873 1875 vx_input->GetParameterValue(&vx,gauss); 1874 1876 vy_input->GetParameterValue(&vy,gauss); … … 1878 1880 1879 1881 /*Loop over all requested responses*/ 1880 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){ 1881 1882 case SurfaceAbsVelMisfitEnum: 1883 /* 1884 * 1 [ 2 2 ] 1885 * J = --- | (u - u ) + (v - v ) | 1886 * 2 [ obs obs ] 1887 * 1888 * dJ 1889 * DU = - -- = (u - u ) 1890 * du obs 1891 */ 1892 for (i=0;i<NUMVERTICES;i++){ 1893 dux=vxobs-vx; 1894 duy=vyobs-vy; 1895 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1896 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1897 } 1898 break; 1899 case SurfaceRelVelMisfitEnum: 1900 /* 1901 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1902 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1903 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1904 * obs obs 1905 * 1906 * dJ \bar{v}^2 1907 * DU = - -- = ------------- (u - u ) 1908 * du (u + eps)^2 obs 1909 * obs 1910 */ 1911 for (i=0;i<NUMVERTICES;i++){ 1912 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1913 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1914 dux=scalex*(vxobs-vx); 1915 duy=scaley*(vyobs-vy); 1916 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1917 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1918 } 1919 break; 1920 case SurfaceLogVelMisfitEnum: 1921 /* 1922 * [ vel + eps ] 2 1923 * J = 4 \bar{v}^2 | log ( ----------- ) | 1924 * [ vel + eps ] 1925 * obs 1926 * 1927 * dJ 2 * log(...) 1928 * DU = - -- = - 4 \bar{v}^2 ------------- u 1929 * du vel^2 + eps 1930 * 1931 */ 1932 for (i=0;i<NUMVERTICES;i++){ 1933 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1934 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1935 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1936 dux=scale*vx; 1937 duy=scale*vy; 1938 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1939 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1940 } 1941 break; 1942 case SurfaceAverageVelMisfitEnum: 1943 /* 1944 * 1 2 2 1945 * J = --- sqrt( (u - u ) + (v - v ) ) 1946 * S obs obs 1947 * 1948 * dJ 1 1 1949 * DU = - -- = - --- ----------- * 2 (u - u ) 1950 * du S 2 sqrt(...) obs 1951 */ 1952 for (i=0;i<NUMVERTICES;i++){ 1953 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1954 dux=scale*(vxobs-vx); 1955 duy=scale*(vyobs-vy); 1956 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1957 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1958 } 1959 break; 1960 case SurfaceLogVxVyMisfitEnum: 1961 /* 1962 * 1 [ |u| + eps 2 |v| + eps 2 ] 1963 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1964 * 2 [ |u |+ eps |v |+ eps ] 1965 * obs obs 1966 * dJ 1 u 1 1967 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1968 * du |u| + eps |u| u + eps 1969 */ 1970 for (i=0;i<NUMVERTICES;i++){ 1971 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1972 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1973 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1974 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1975 } 1976 break; 1977 default: 1978 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1882 for(resp=0;resp<num_responses;resp++){ 1883 1884 weights_input->GetParameterValue(&weight,gauss,resp); 1885 1886 switch(responses[resp]){ 1887 1888 case SurfaceAbsVelMisfitEnum: 1889 /* 1890 * 1 [ 2 2 ] 1891 * J = --- | (u - u ) + (v - v ) | 1892 * 2 [ obs obs ] 1893 * 1894 * dJ 1895 * DU = - -- = (u - u ) 1896 * du obs 1897 */ 1898 for (i=0;i<NUMVERTICES;i++){ 1899 dux=vxobs-vx; 1900 duy=vyobs-vy; 1901 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1902 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1903 } 1904 break; 1905 case SurfaceRelVelMisfitEnum: 1906 /* 1907 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1908 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1909 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1910 * obs obs 1911 * 1912 * dJ \bar{v}^2 1913 * DU = - -- = ------------- (u - u ) 1914 * du (u + eps)^2 obs 1915 * obs 1916 */ 1917 for (i=0;i<NUMVERTICES;i++){ 1918 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0; 1919 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0; 1920 dux=scalex*(vxobs-vx); 1921 duy=scaley*(vyobs-vy); 1922 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1923 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1924 } 1925 break; 1926 case SurfaceLogVelMisfitEnum: 1927 /* 1928 * [ vel + eps ] 2 1929 * J = 4 \bar{v}^2 | log ( ----------- ) | 1930 * [ vel + eps ] 1931 * obs 1932 * 1933 * dJ 2 * log(...) 1934 * DU = - -- = - 4 \bar{v}^2 ------------- u 1935 * du vel^2 + eps 1936 * 1937 */ 1938 for (i=0;i<NUMVERTICES;i++){ 1939 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel; 1940 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel; 1941 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag); 1942 dux=scale*vx; 1943 duy=scale*vy; 1944 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1945 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1946 } 1947 break; 1948 case SurfaceAverageVelMisfitEnum: 1949 /* 1950 * 1 2 2 1951 * J = --- sqrt( (u - u ) + (v - v ) ) 1952 * S obs obs 1953 * 1954 * dJ 1 1 1955 * DU = - -- = - --- ----------- * 2 (u - u ) 1956 * du S 2 sqrt(...) obs 1957 */ 1958 for (i=0;i<NUMVERTICES;i++){ 1959 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel); 1960 dux=scale*(vxobs-vx); 1961 duy=scale*(vyobs-vy); 1962 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1963 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1964 } 1965 break; 1966 case SurfaceLogVxVyMisfitEnum: 1967 /* 1968 * 1 [ |u| + eps 2 |v| + eps 2 ] 1969 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1970 * 2 [ |u |+ eps |v |+ eps ] 1971 * obs obs 1972 * dJ 1 u 1 1973 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1974 * du |u| + eps |u| u + eps 1975 */ 1976 for (i=0;i<NUMVERTICES;i++){ 1977 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1978 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1979 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 1980 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 1981 } 1982 break; 1983 default: 1984 _error_("response %s not supported yet",EnumToStringx(responses[resp])); 1985 } 1979 1986 } 1980 1987 } … … 4820 4827 /*}}}*/ 4821 4828 /*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ 4822 double Tria::SurfaceAbsVelMisfit(bool process_units ){4829 double Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){ 4823 4830 4824 4831 const int numdof=NDOF2*NUMVERTICES; … … 4854 4861 4855 4862 /*Get all parameters at gaussian point*/ 4856 weights_input->GetParameterValue(&weight,gauss, 0);4863 weights_input->GetParameterValue(&weight,gauss,weight_index); 4857 4864 vx_input->GetParameterValue(&vx,gauss); 4858 4865 vy_input->GetParameterValue(&vy,gauss); … … 4910 4917 /*}}}*/ 4911 4918 /*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/ 4912 double Tria::SurfaceAverageVelMisfit(bool process_units ){4919 double Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){ 4913 4920 4914 4921 const int numdof=2*NUMVERTICES; … … 4945 4952 4946 4953 /*Get all parameters at gaussian point*/ 4947 weights_input->GetParameterValue(&weight,gauss, 0);4954 weights_input->GetParameterValue(&weight,gauss,weight_index); 4948 4955 vx_input->GetParameterValue(&vx,gauss); 4949 4956 vy_input->GetParameterValue(&vy,gauss); … … 4971 4978 /*}}}*/ 4972 4979 /*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/ 4973 double Tria::SurfaceLogVelMisfit(bool process_units ){4980 double Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){ 4974 4981 4975 4982 const int numdof=NDOF2*NUMVERTICES; … … 5008 5015 5009 5016 /*Get all parameters at gaussian point*/ 5010 weights_input->GetParameterValue(&weight,gauss, 0);5017 weights_input->GetParameterValue(&weight,gauss,weight_index); 5011 5018 vx_input->GetParameterValue(&vx,gauss); 5012 5019 vy_input->GetParameterValue(&vy,gauss); … … 5036 5043 /*}}}*/ 5037 5044 /*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/ 5038 double Tria::SurfaceLogVxVyMisfit(bool process_units ){5045 double Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){ 5039 5046 5040 5047 const int numdof=NDOF2*NUMVERTICES; … … 5073 5080 5074 5081 /*Get all parameters at gaussian point*/ 5075 weights_input->GetParameterValue(&weight,gauss, 0);5082 weights_input->GetParameterValue(&weight,gauss,weight_index); 5076 5083 vx_input->GetParameterValue(&vx,gauss); 5077 5084 vy_input->GetParameterValue(&vy,gauss); … … 5126 5133 /*}}}*/ 5127 5134 /*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ 5128 double Tria::SurfaceRelVelMisfit(bool process_units ){5135 double Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){ 5129 5136 const int numdof=2*NUMVERTICES; 5130 5137 … … 5162 5169 5163 5170 /*Get all parameters at gaussian point*/ 5164 weights_input->GetParameterValue(&weight,gauss, 0);5171 weights_input->GetParameterValue(&weight,gauss,weight_index); 5165 5172 vx_input->GetParameterValue(&vx,gauss); 5166 5173 vy_input->GetParameterValue(&vy,gauss); … … 5190 5197 /*}}}*/ 5191 5198 /*FUNCTION Tria::ThicknessAbsMisfit {{{1*/ 5192 double Tria::ThicknessAbsMisfit(bool process_units ){5199 double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){ 5193 5200 5194 5201 /* Constants */ … … 5230 5237 thickness_input->GetParameterDerivativeValue(&dH[0],&xyz_list[0][0],gauss); 5231 5238 thicknessobs_input->GetParameterValue(&thicknessobs,gauss); 5232 weights_input->GetParameterValue(&weight,gauss, 0);5239 weights_input->GetParameterValue(&weight,gauss,weight_index); 5233 5240 5234 5241 /*compute ThicknessAbsMisfit*/ -
issm/trunk/src/c/objects/Elements/Tria.h
r8592 r8607 123 123 void MinVy(double* pminvy, bool process_units); 124 124 void MinVz(double* pminvz, bool process_units); 125 double ThicknessAbsMisfit( bool process_units);126 double SurfaceAbsVelMisfit( bool process_units);127 double SurfaceRelVelMisfit( bool process_units);128 double SurfaceLogVelMisfit( bool process_units);129 double SurfaceLogVxVyMisfit( bool process_units);130 double SurfaceAverageVelMisfit(bool process_units );125 double ThicknessAbsMisfit( bool process_units,int weight_index); 126 double SurfaceAbsVelMisfit( bool process_units,int weight_index); 127 double SurfaceRelVelMisfit( bool process_units,int weight_index); 128 double SurfaceLogVelMisfit( bool process_units,int weight_index); 129 double SurfaceLogVxVyMisfit( bool process_units,int weight_index); 130 double SurfaceAverageVelMisfit(bool process_units,int weight_index); 131 131 void PatchFill(int* pcount, Patch* patch); 132 132 void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); -
issm/trunk/src/c/solutions/controltao_core.cpp
r8429 r8607 113 113 //VecScale(gradient,-1.); 114 114 VecCopy(gradient,G); 115 CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters ,SurfaceAbsVelMisfitEnum);115 CostFunctionx(fcn, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 116 116 117 117 //printf("X\n"); -
issm/trunk/src/c/solutions/objectivefunctionC.cpp
r8601 r8607 25 25 26 26 /*output: */ 27 double J =0,Jplus;27 double J; 28 28 29 29 /*parameters: */ 30 int num_responses;31 30 int solution_type,analysis_type; 32 31 bool isstokes = false; 33 32 bool conserve_loads = true; 34 int *responses = NULL;35 33 FemModel *femmodel = NULL; 36 34 … … 39 37 40 38 /*Recover parameters: */ 41 femmodel->parameters->FindParam(&num_responses,NumResponsesEnum);42 femmodel->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);43 39 femmodel->parameters->FindParam(&isstokes,IsStokesEnum); 44 40 femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 74 70 75 71 /*Compute misfit for this velocity field.*/ 76 for(int i=0;i<num_responses;i++){ 77 CostFunctionx(&Jplus, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,responses[i]); 78 J+=Jplus; 79 } 72 CostFunctionx(&J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 80 73 81 74 /*Free ressources:*/ 82 xfree((void**)&responses);83 75 return J; 84 76 } -
issm/trunk/src/m/solutions/objectivefunctionC.m
r8602 r8607 6 6 7 7 %recover some parameters 8 num_responses = femmodel.parameters.NumResponses;9 responses = femmodel.parameters.StepResponses;10 8 analysis_type = femmodel.parameters.AnalysisType; 11 9 solution_type = femmodel.parameters.SolutionType; … … 36 34 37 35 %Compute misfit for this velocity field 38 for i=1:num_responses 39 J=J+CostFunction(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials, femmodel.parameters,responses(i)); 40 end 36 J=CostFunction(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials, femmodel.parameters); -
issm/trunk/src/mex/CostFunction/CostFunction.cpp
r6412 r8607 14 14 Materials *materials = NULL; 15 15 Parameters *parameters = NULL; 16 int response;17 16 18 17 /* output datasets: */ … … 32 31 FetchData((DataSet**)&materials,MATERIALS); 33 32 FetchParams(¶meters,PARAMETERS); 34 FetchData(&response,RESPONSE);35 33 36 34 /*configure: */ … … 40 38 41 39 /*!Call core code: */ 42 CostFunctionx(&J, elements,nodes,vertices, loads,materials,parameters ,response);40 CostFunctionx(&J, elements,nodes,vertices, loads,materials,parameters); 43 41 44 42 /*write output : */ … … 60 58 { 61 59 _printf_(true,"\n"); 62 _printf_(true," usage: [J] = %s(elements,nodes,vertices,loads, materials, parameters ,response);\n",__FUNCT__);60 _printf_(true," usage: [J] = %s(elements,nodes,vertices,loads, materials, parameters);\n",__FUNCT__); 63 61 _printf_(true,"\n"); 64 62 } -
issm/trunk/src/mex/CostFunction/CostFunction.h
r5280 r8607 23 23 #define MATERIALS (mxArray*)prhs[4] 24 24 #define PARAMETERS (mxArray*)prhs[5] 25 #define RESPONSE (mxArray*)prhs[6]26 25 27 26 /* serial output macros: */ … … 32 31 #define NLHS 1 33 32 #undef NRHS 34 #define NRHS 733 #define NRHS 6 35 34 36 35 #endif /* _COSTFUNCTION_H */ -
issm/trunk/src/mex/Response/Response.cpp
r6716 r8607 16 16 char *response = NULL; 17 17 bool process_units; 18 int weight_index; 18 19 19 20 /* output datasets: */ … … 35 36 FetchData(&response,RESPONSE); 36 37 FetchData(&process_units,PROCESSUNITS); 38 FetchData(&weight_index,WEIGHTINDEX); 37 39 38 40 /*configure: */ … … 42 44 43 45 /*!Call core code: */ 44 Responsex(&resp, elements,nodes,vertices, loads,materials,parameters,response,process_units );46 Responsex(&resp, elements,nodes,vertices, loads,materials,parameters,response,process_units,weight_index); 45 47 46 48 /*write output : */ -
issm/trunk/src/mex/Response/Response.h
r6689 r8607 26 26 #define RESPONSE (mxArray*)prhs[6] 27 27 #define PROCESSUNITS (mxArray*)prhs[7] 28 #define WEIGHTINDEX (mxArray*)prhs[8] 28 29 29 30 /* serial output macros: */ … … 34 35 #define NLHS 1 35 36 #undef NRHS 36 #define NRHS 837 #define NRHS 9 37 38 38 39 #endif /* _RESPONSE_H */
Note:
See TracChangeset
for help on using the changeset viewer.