![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5281
- Timestamp:
- Aug 16, 2010, 11:37:04 AM (16 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 1 deleted
- 24 edited
-
EnumDefinitions/EnumDefinitions.h (modified) (3 diffs)
-
EnumDefinitions/EnumToString.cpp (modified) (3 diffs)
-
EnumDefinitions/StringToEnum.cpp (modified) (4 diffs)
-
Makefile.am (modified) (2 diffs)
-
modules/CostFunctionx/CostFunctionx.cpp (modified) (2 diffs)
-
modules/CostFunctionx/CostFunctionx.h (modified) (1 diff)
-
modules/MassFluxx/MassFluxx.cpp (modified) (1 diff)
-
modules/MassFluxx/MassFluxx.h (modified) (1 diff)
-
modules/Misfitx (deleted)
-
modules/ModelProcessorx/Control/CreateParametersControl.cpp (modified) (1 diff)
-
modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp (modified) (1 diff)
-
modules/Responsex/Responsex.cpp (modified) (2 diffs)
-
modules/Responsex/Responsex.h (modified) (1 diff)
-
modules/modules.h (modified) (1 diff)
-
objects/Elements/Element.h (modified) (1 diff)
-
objects/Elements/Penta.cpp (modified) (8 diffs)
-
objects/Elements/Penta.h (modified) (2 diffs)
-
objects/Elements/Tria.cpp (modified) (27 diffs)
-
objects/Elements/Tria.h (modified) (2 diffs)
-
objects/IoModel.cpp (modified) (2 diffs)
-
objects/IoModel.h (modified) (1 diff)
-
solutions/control_core.cpp (modified) (6 diffs)
-
solutions/controlconvergence.cpp (modified) (2 diffs)
-
solutions/objectivefunctionC.cpp (modified) (3 diffs)
-
solutions/solutions.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
r5251 r5281 170 170 BedSlopeYEnum, 171 171 BoundaryEnum, 172 CmResponseEnum, 173 CmResponsesEnum, 172 174 CmMaxDmpSlopeEnum, 173 175 CmMaxDmpValueEnum, … … 190 192 EpsVelEnum, 191 193 FillEnum, 192 FitEnum,193 194 FractionIncrementEnum, 194 195 FrictionEnum, … … 197 198 KflagEnum, 198 199 MassFluxEnum, 200 SurfaceAbsVelMisfitEnum, 201 SurfaceRelVelMisfitEnum, 202 SurfaceLogVelMisfitEnum, 203 SurfaceLogVxVyMisfitEnum, 204 SurfaceAverageVelMisfitEnum, 199 205 MaxPenetrationEnum, 200 206 MeanVelEnum, -
issm/trunk/src/c/EnumDefinitions/EnumToString.cpp
r5251 r5281 144 144 case BedSlopeYEnum : return "BedSlopeY"; 145 145 case BoundaryEnum : return "Boundary"; 146 case CmResponseEnum : return "CmResponse"; 147 case CmResponsesEnum : return "CmResponses"; 146 148 case CmMaxDmpSlopeEnum : return "CmMaxDmpSlope"; 147 149 case CmMaxDmpValueEnum : return "CmMaxDmpValue"; … … 164 166 case EpsVelEnum : return "EpsVel"; 165 167 case FillEnum : return "Fill"; 166 case FitEnum : return "Fit";167 168 case FractionIncrementEnum : return "FractionIncrement"; 168 169 case FrictionEnum : return "Friction"; … … 171 172 case KflagEnum : return "Kflag"; 172 173 case MassFluxEnum : return "MassFlux"; 174 case SurfaceAbsVelMisfitEnum : return "SurfaceAbsVelMisfit"; 175 case SurfaceRelVelMisfitEnum : return "SurfaceRelVelMisfit"; 176 case SurfaceLogVelMisfitEnum : return "SurfaceLogVelMisfit"; 177 case SurfaceLogVxVyMisfitEnum : return "SurfaceLogVxVyMisfit"; 178 case SurfaceAverageVelMisfitEnum : return "SurfaceAverageVelMisfit"; 173 179 case MaxPenetrationEnum : return "MaxPenetration"; 174 180 case MeanVelEnum : return "MeanVel"; -
issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp
r5251 r5281 142 142 else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; 143 143 else if (strcmp(name,"Boundary")==0) return BoundaryEnum; 144 else if (strcmp(name,"CmResponse")==0) return CmResponseEnum; 145 else if (strcmp(name,"CmResponses")==0) return CmResponsesEnum; 144 146 else if (strcmp(name,"CmMaxDmpSlope")==0) return CmMaxDmpSlopeEnum; 145 147 else if (strcmp(name,"CmMaxDmpValue")==0) return CmMaxDmpValueEnum; … … 162 164 else if (strcmp(name,"EpsVel")==0) return EpsVelEnum; 163 165 else if (strcmp(name,"Fill")==0) return FillEnum; 164 else if (strcmp(name,"Fit")==0) return FitEnum;165 166 else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum; 166 167 else if (strcmp(name,"Friction")==0) return FrictionEnum; … … 169 170 else if (strcmp(name,"Kflag")==0) return KflagEnum; 170 171 else if (strcmp(name,"MassFlux")==0) return MassFluxEnum; 172 else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; 173 else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; 174 else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; 175 else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; 176 else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; 171 177 else if (strcmp(name,"MaxPenetration")==0) return MaxPenetrationEnum; 172 178 else if (strcmp(name,"MeanVel")==0) return MeanVelEnum; … … 305 311 else if (strcmp(name,"WaitOnLock")==0) return WaitOnLockEnum; 306 312 else if (strcmp(name,"Yts")==0) return YtsEnum; 307 else ISSMERROR("Enum %s not found",name);308 309 } -
issm/trunk/src/c/Makefile.am
r5250 r5281 433 433 ./modules/SurfaceAreax/SurfaceAreax.h\ 434 434 ./modules/SurfaceAreax/SurfaceAreax.cpp\ 435 ./modules/Misfitx/Misfitx.h\ 436 ./modules/Misfitx/Misfitx.cpp\ 435 ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h\ 436 ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp\ 437 ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h\ 438 ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp\ 439 ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h\ 440 ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp\ 441 ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h\ 442 ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp\ 443 ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\ 444 ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\ 437 445 ./modules/CostFunctionx/CostFunctionx.h\ 438 446 ./modules/CostFunctionx/CostFunctionx.cpp\ … … 968 976 ./modules/SurfaceAreax/SurfaceAreax.h\ 969 977 ./modules/SurfaceAreax/SurfaceAreax.cpp\ 970 ./modules/Misfitx/Misfitx.h\ 971 ./modules/Misfitx/Misfitx.cpp\ 978 ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h\ 979 ./modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp\ 980 ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h\ 981 ./modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp\ 982 ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h\ 983 ./modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp\ 984 ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h\ 985 ./modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.cpp\ 986 ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\ 987 ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\ 972 988 ./modules/CostFunctionx/CostFunctionx.h\ 973 989 ./modules/CostFunctionx/CostFunctionx.cpp\ -
issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.cpp
r4974 r5281 9 9 #include "../../toolkits/toolkits.h" 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 #include "../ SurfaceAreax/SurfaceAreax.h"11 #include "../Responsex/Responsex.h" 12 12 13 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters, bool process_units){13 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters,int response){ 14 14 15 15 /*Intermediary*/ … … 20 20 21 21 /*output: */ 22 double J=0; 23 double J_sum; 22 double J; 23 double Jreg=0; 24 double Jreg_sum; 25 26 /*Get response*/ 27 Responsex(&J,elements,nodes,vertices,loads,materials,parameters,(const char*)EnumToString(response),false); //False means DO NOT process units 24 28 25 /*Compute surface area: */ 26 SurfaceAreax(&S,elements,nodes,vertices, loads,materials,parameters); 27 28 /*add surface area to elements :*/ 29 for(i=0;i<elements->Size();i++){ 30 Element* element=(Element*)elements->GetObjectByOffset(i); 31 element->InputUpdateFromVector(&S,SurfaceAreaEnum,ConstantEnum); 32 } 33 34 /*Compute cost function: */ 29 /*Add Regularization terms: */ 35 30 for (i=0;i<elements->Size();i++){ 36 31 element=(Element*)elements->GetObjectByOffset(i); 37 J +=element->CostFunction(process_units);32 Jreg+=element->RegularizeInversion(); 38 33 } 39 34 40 35 /*Sum all J from all cpus of the cluster:*/ 41 MPI_Reduce (&J ,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD );42 MPI_Bcast(&J _sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD);43 J=J _sum;36 MPI_Reduce (&Jreg,&Jreg_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD ); 37 MPI_Bcast(&Jreg_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 38 J=J+Jreg_sum; 44 39 45 40 /*Assign output pointers: */ -
issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h
r4974 r5281 10 10 11 11 /* local prototypes: */ 12 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters, bool process_units);12 void CostFunctionx( double* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,int response); 13 13 14 #endif /* _MISFITX_H */ 15 14 #endif -
issm/trunk/src/c/modules/MassFluxx/MassFluxx.cpp
r5250 r5281 10 10 #include "../../EnumDefinitions/EnumDefinitions.h" 11 11 12 void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,c har* descriptor,bool process_units){12 void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* descriptor,bool process_units){ 13 13 14 14 int i,j; -
issm/trunk/src/c/modules/MassFluxx/MassFluxx.h
r5250 r5281 10 10 11 11 /* local prototypes: */ 12 void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,c har* descriptor,bool process_units);12 void MassFluxx(double* pmass_flux, Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* descriptor,bool process_units); 13 13 14 14 -
issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
r5247 r5281 60 60 61 61 /*Now, recover fit, optscal and maxiter as vectors: */ 62 IoModelFetchData(&iomodel-> fit,NULL,NULL,iomodel_handle,"fit");62 IoModelFetchData(&iomodel->cm_responses,NULL,NULL,iomodel_handle,"cm_responses"); 63 63 IoModelFetchData(&iomodel->cm_jump,NULL,NULL,iomodel_handle,"cm_jump"); 64 64 IoModelFetchData(&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal"); 65 65 IoModelFetchData(&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter"); 66 66 67 parameters->AddObject(new DoubleVecParam( FitEnum,iomodel->fit,iomodel->nsteps));67 parameters->AddObject(new DoubleVecParam(CmResponsesEnum,iomodel->cm_responses,iomodel->nsteps)); 68 68 parameters->AddObject(new DoubleVecParam(CmJumpEnum,iomodel->cm_jump,iomodel->nsteps)); 69 69 parameters->AddObject(new DoubleVecParam(OptScalEnum,iomodel->optscal,iomodel->nsteps)); 70 70 parameters->AddObject(new DoubleVecParam(MaxIterEnum,iomodel->maxiter,iomodel->nsteps)); 71 71 72 xfree((void**)&iomodel-> fit);72 xfree((void**)&iomodel->cm_responses); 73 73 xfree((void**)&iomodel->cm_jump); 74 74 xfree((void**)&iomodel->optscal); -
issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp
r5240 r5281 207 207 qmu_mass_flux_present=true; 208 208 } 209 210 if (strcmp(descriptor,"misfit")==0){211 212 /*We need the fit: */213 IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit");214 parameters->SetParam(iomodel->fit,iomodel->nsteps,FitEnum);215 xfree((void**)&iomodel->fit);216 217 }218 209 } 219 210 -
issm/trunk/src/c/modules/Responsex/Responsex.cpp
r5250 r5281 16 16 #include "../modules.h" 17 17 18 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,c har* response_descriptor,bool process_units){18 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* response_descriptor,bool process_units){ 19 19 20 20 if(strcmp(response_descriptor,"MinVel")==0){ … … 51 51 MaxAbsVzx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); 52 52 } 53 else if(strcmp(response_descriptor," Misfit")==0){54 Misfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);53 else if(strcmp(response_descriptor,"SurfaceAbsVelMisfit")==0){ 54 SurfaceAbsVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); 55 55 } 56 else if(strlen(response_descriptor)>=8){ 57 if(strncmp(response_descriptor,"MassFlux",8)==0) MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units); 56 else if(strcmp(response_descriptor,"SurfaceRelVelMisfit")==0){ 57 SurfaceRelVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); 58 } 59 else if(strcmp(response_descriptor,"SurfaceLogVelMisfit")==0){ 60 SurfaceLogVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); 61 } 62 else if(strcmp(response_descriptor,"SurfaceLogVxVyMisfit")==0){ 63 SurfaceLogVxVyMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); 64 } 65 else if(strcmp(response_descriptor,"SurfaceAverageVelMisfit")==0){ 66 SurfaceAverageVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units); 67 } 68 else if(strncmp(response_descriptor,"MassFlux",8)==0){ 69 MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units); 58 70 } 59 71 else{ 60 ISSMERROR(" %s%s%s"," response descriptor : ",response_descriptor," not supported yet!");72 ISSMERROR(" response descriptor \"%s\" not supported yet!",response_descriptor); 61 73 } 62 74 } -
issm/trunk/src/c/modules/Responsex/Responsex.h
r5250 r5281 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,c har* response_descriptor,bool process_units);11 void Responsex(double* presponse,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,const char* response_descriptor,bool process_units); 12 12 13 13 #endif /* _RESPONSESXX_H */ -
issm/trunk/src/c/modules/modules.h
r5250 r5281 55 55 #include "./MinVyx/MinVyx.h" 56 56 #include "./MinVzx/MinVzx.h" 57 #include "./Misfitx/Misfitx.h" 57 #include "./SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h" 58 #include "./SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h" 59 #include "./SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h" 60 #include "./SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h" 61 #include "./SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h" 58 62 #include "./ModelProcessorx/ModelProcessorx.h" 59 63 #include "./NodeConnectivityx/NodeConnectivityx.h" -
issm/trunk/src/c/objects/Elements/Element.h
r4974 r5281 42 42 virtual void GradjDrag(Vec gradient)=0; 43 43 virtual void GradjB(Vec gradient)=0; 44 virtual double Misfit(bool process_units)=0; 45 virtual double CostFunction(bool process_units)=0; 44 virtual double SurfaceAbsVelMisfit(bool process_units)=0; 45 virtual double SurfaceRelVelMisfit(bool process_units)=0; 46 virtual double SurfaceLogVelMisfit(bool process_units)=0; 47 virtual double SurfaceLogVxVyMisfit(bool process_units)=0; 48 virtual double SurfaceAverageVelMisfit(bool process_units)=0; 49 virtual double RegularizeInversion(void)=0; 46 50 virtual double SurfaceArea(void)=0; 47 51 virtual void InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0; -
issm/trunk/src/c/objects/Elements/Penta.cpp
r5263 r5281 666 666 } 667 667 /*}}}*/ 668 /*FUNCTION Penta:: CostFunction {{{1*/669 double Penta:: CostFunction(bool process_units){668 /*FUNCTION Penta::RegularizeInversion {{{1*/ 669 double Penta::RegularizeInversion(void){ 670 670 671 671 double J; … … 696 696 697 697 /*This element should be collapsed into a tria element at its base. Create this tria element, 698 * and compute CostFunction*/698 * and compute RegularizeInversion*/ 699 699 700 700 /*Depth Average B*/ … … 702 702 703 703 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 704 J=tria-> CostFunction(process_units);704 J=tria->RegularizeInversion(); 705 705 delete tria->matice; delete tria; 706 706 … … 720 720 721 721 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 722 J=tria-> CostFunction(process_units);722 J=tria->RegularizeInversion(); 723 723 delete tria->matice; delete tria; 724 724 … … 1597 1597 } 1598 1598 /*}}}*/ 1599 /*FUNCTION Penta:: Misfit {{{1*/1600 double Penta:: Misfit(bool process_units){1599 /*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/ 1600 double Penta::SurfaceAbsVelMisfit(bool process_units){ 1601 1601 1602 1602 double J; … … 1627 1627 1628 1628 /*This element should be collapsed into a tria element at its base. Create this tria element, 1629 * and compute Misfit*/1629 * and compute SurfaceAbsVelMisfit*/ 1630 1630 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 1631 J=tria-> Misfit(process_units);1631 J=tria->SurfaceAbsVelMisfit(process_units); 1632 1632 delete tria->matice; delete tria; 1633 1633 return J; … … 1636 1636 1637 1637 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 1638 J=tria->Misfit(process_units); 1638 J=tria->SurfaceAbsVelMisfit(process_units); 1639 delete tria->matice; delete tria; 1640 return J; 1641 } 1642 } 1643 /*}}}*/ 1644 /*FUNCTION Penta::SurfaceRelVelMisfit {{{1*/ 1645 double Penta::SurfaceRelVelMisfit(bool process_units){ 1646 1647 double J; 1648 Tria* tria=NULL; 1649 1650 /*inputs: */ 1651 bool onwater; 1652 bool onsurface; 1653 bool onbed; 1654 int approximation; 1655 1656 /*retrieve inputs :*/ 1657 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1658 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1659 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 1660 inputs->GetParameterValue(&approximation,ApproximationEnum); 1661 1662 /*If on water, return 0: */ 1663 if(onwater)return 0; 1664 1665 /*Bail out if this element if: 1666 * -> Non MacAyeal and not on the surface 1667 * -> MacAyeal (2d model) and not on bed) */ 1668 if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ 1669 return 0; 1670 } 1671 else if (MacAyealApproximationEnum){ 1672 1673 /*This element should be collapsed into a tria element at its base. Create this tria element, 1674 * and compute SurfaceRelVelMisfit*/ 1675 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 1676 J=tria->SurfaceRelVelMisfit(process_units); 1677 delete tria->matice; delete tria; 1678 return J; 1679 } 1680 else{ 1681 1682 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 1683 J=tria->SurfaceRelVelMisfit(process_units); 1684 delete tria->matice; delete tria; 1685 return J; 1686 } 1687 } 1688 /*}}}*/ 1689 /*FUNCTION Penta::SurfaceLogVelMisfit {{{1*/ 1690 double Penta::SurfaceLogVelMisfit(bool process_units){ 1691 1692 double J; 1693 Tria* tria=NULL; 1694 1695 /*inputs: */ 1696 bool onwater; 1697 bool onsurface; 1698 bool onbed; 1699 int approximation; 1700 1701 /*retrieve inputs :*/ 1702 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1703 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1704 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 1705 inputs->GetParameterValue(&approximation,ApproximationEnum); 1706 1707 /*If on water, return 0: */ 1708 if(onwater)return 0; 1709 1710 /*Bail out if this element if: 1711 * -> Non MacAyeal and not on the surface 1712 * -> MacAyeal (2d model) and not on bed) */ 1713 if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ 1714 return 0; 1715 } 1716 else if (MacAyealApproximationEnum){ 1717 1718 /*This element should be collapsed into a tria element at its base. Create this tria element, 1719 * and compute SurfaceLogVelMisfit*/ 1720 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 1721 J=tria->SurfaceLogVelMisfit(process_units); 1722 delete tria->matice; delete tria; 1723 return J; 1724 } 1725 else{ 1726 1727 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 1728 J=tria->SurfaceLogVelMisfit(process_units); 1729 delete tria->matice; delete tria; 1730 return J; 1731 } 1732 } 1733 /*}}}*/ 1734 /*FUNCTION Penta::SurfaceLogVxVyMisfit {{{1*/ 1735 double Penta::SurfaceLogVxVyMisfit(bool process_units){ 1736 1737 double J; 1738 Tria* tria=NULL; 1739 1740 /*inputs: */ 1741 bool onwater; 1742 bool onsurface; 1743 bool onbed; 1744 int approximation; 1745 1746 /*retrieve inputs :*/ 1747 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1748 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1749 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 1750 inputs->GetParameterValue(&approximation,ApproximationEnum); 1751 1752 /*If on water, return 0: */ 1753 if(onwater)return 0; 1754 1755 /*Bail out if this element if: 1756 * -> Non MacAyeal and not on the surface 1757 * -> MacAyeal (2d model) and not on bed) */ 1758 if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ 1759 return 0; 1760 } 1761 else if (MacAyealApproximationEnum){ 1762 1763 /*This element should be collapsed into a tria element at its base. Create this tria element, 1764 * and compute SurfaceLogVxVyMisfit*/ 1765 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 1766 J=tria->SurfaceLogVxVyMisfit(process_units); 1767 delete tria->matice; delete tria; 1768 return J; 1769 } 1770 else{ 1771 1772 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 1773 J=tria->SurfaceLogVxVyMisfit(process_units); 1774 delete tria->matice; delete tria; 1775 return J; 1776 } 1777 } 1778 /*}}}*/ 1779 /*FUNCTION Penta::SurfaceAverageVelMisfit {{{1*/ 1780 double Penta::SurfaceAverageVelMisfit(bool process_units){ 1781 1782 double J; 1783 Tria* tria=NULL; 1784 1785 /*inputs: */ 1786 bool onwater; 1787 bool onsurface; 1788 bool onbed; 1789 int approximation; 1790 1791 /*retrieve inputs :*/ 1792 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 1793 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 1794 inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); 1795 inputs->GetParameterValue(&approximation,ApproximationEnum); 1796 1797 /*If on water, return 0: */ 1798 if(onwater)return 0; 1799 1800 /*Bail out if this element if: 1801 * -> Non MacAyeal and not on the surface 1802 * -> MacAyeal (2d model) and not on bed) */ 1803 if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){ 1804 return 0; 1805 } 1806 else if (MacAyealApproximationEnum){ 1807 1808 /*This element should be collapsed into a tria element at its base. Create this tria element, 1809 * and compute SurfaceAverageVelMisfit*/ 1810 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). 1811 J=tria->SurfaceAverageVelMisfit(process_units); 1812 delete tria->matice; delete tria; 1813 return J; 1814 } 1815 else{ 1816 1817 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). 1818 J=tria->SurfaceAverageVelMisfit(process_units); 1639 1819 delete tria->matice; delete tria; 1640 1820 return J; … … 5900 6080 name==TemperatureEnum || 5901 6081 name==ControlParameterEnum || 5902 name== FitEnum ||6082 name==CmResponseEnum || 5903 6083 name==DragCoefficientEnum || 5904 6084 name==GradientEnum || -
issm/trunk/src/c/objects/Elements/Penta.h
r5258 r5281 71 71 void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); 72 72 void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); 73 double CostFunction(bool process_units);73 double RegularizeInversion(void); 74 74 void CreateKMatrix(Mat Kgg); 75 75 void CreatePVector(Vec pg); … … 104 104 void MinVy(double* pminvy, bool process_units); 105 105 void MinVz(double* pminvz, bool process_units); 106 double Misfit(bool process_units); 106 double SurfaceAbsVelMisfit(bool process_units); 107 double SurfaceRelVelMisfit(bool process_units); 108 double SurfaceLogVelMisfit(bool process_units); 109 double SurfaceLogVxVyMisfit(bool process_units); 110 double SurfaceAverageVelMisfit(bool process_units); 107 111 void PatchFill(int* pcount, Patch* patch); 108 112 void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r5263 r5281 548 548 } 549 549 /*}}}*/ 550 /*FUNCTION Tria:: CostFunction {{{1*/551 double Tria:: CostFunction(bool process_units){550 /*FUNCTION Tria::RegularizeInversion {{{1*/ 551 double Tria::RegularizeInversion(){ 552 552 553 553 int i; 554 554 555 555 /* output: */ 556 double Jelem ;556 double Jelem=0; 557 557 558 558 /* node data: */ … … 602 602 /* Get node coordinates and dof list: */ 603 603 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 604 605 /*First, get Misfit*/606 Jelem=Misfit(process_units);607 604 608 605 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 1860 1857 } 1861 1858 /*}}}*/ 1862 /*FUNCTION Tria:: Misfit {{{1*/1863 double Tria:: Misfit(bool process_units){1859 /*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ 1860 double Tria::SurfaceAbsVelMisfit(bool process_units){ 1864 1861 1865 1862 int i; … … 1900 1897 double Jdet; 1901 1898 1902 /*relative and logarithmic control method :*/1903 double scalex=1;1904 double scaley=1;1905 double S=0;1906 int fit=-1;1907 1908 1899 /*inputs: */ 1909 1900 bool onwater; … … 1921 1912 1922 1913 /* Recover input data: */ 1923 inputs->GetParameterValue(&fit,FitEnum);1924 if(fit==3){1925 inputs->GetParameterValue(&S,SurfaceAreaEnum);1926 }1927 1914 inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); 1928 1915 inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); … … 1935 1922 this->parameters->FindParam(&epsvel,EpsVelEnum); 1936 1923 1937 /* Compute Misfit at the 3 nodes1924 /* Compute SurfaceAbsVelMisfit at the 3 nodes 1938 1925 * Here we integrate linearized functions: 1939 1926 * … … 1944 1931 * the vertex i 1945 1932 */ 1946 if(fit==0){ 1947 /*We are using an absolute misfit: 1948 * 1949 * 1 [ 2 2 ] 1950 * J = --- | (u - u ) + (v - v ) | 1951 * 2 [ obs obs ] 1952 * 1953 */ 1954 for (i=0;i<numgrids;i++){ 1955 misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2)); 1956 } 1957 /*Process units: */ 1958 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters); 1959 1960 } 1961 else if(fit==1){ 1962 /*We are using a relative misfit: 1963 * 1964 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1965 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1966 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1967 * obs obs 1968 */ 1969 for (i=0;i<numgrids;i++){ 1970 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2); 1971 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2); 1972 if(obs_vx_list[i]==0)scalex=0; 1973 if(obs_vy_list[i]==0)scaley=0; 1974 misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2)); 1975 } 1976 1977 /*Process units: */ 1978 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters); 1979 1980 } 1981 else if(fit==2){ 1982 /*We are using a logarithmic misfit: 1983 * 1984 * [ vel + eps ] 2 1985 * J = 4 \bar{v}^2 | log ( ----------- ) | 1986 * [ vel + eps ] 1987 * obs 1988 */ 1989 for (i=0;i<numgrids;i++){ 1990 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil. 1991 obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil. 1992 misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); 1993 } 1994 1995 /*Process units: */ 1996 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters); 1997 1998 } 1999 else if(fit==3){ 2000 /*We are using a spacially average absolute misfit: 2001 * 2002 * 1 2 2 2003 * J = --- sqrt( (u - u ) + (v - v ) ) 2004 * S obs obs 2005 */ 2006 for (i=0;i<numgrids;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2); 2007 2008 /*Process units: */ 2009 if(process_units)NodalValuesUnitConversion(&misfit_square_list[0],numgrids,MisfitEnum,this->parameters); 2010 2011 /*Take the square root, and scale by surface: */ 2012 for (i=0;i<numgrids;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S; 2013 } 2014 else if(fit==4){ 2015 /*We are using an logarithmic 2 misfit: 2016 * 2017 * 1 [ |u| + eps 2 |v| + eps 2 ] 2018 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 2019 * 2 [ |u |+ eps |v |+ eps ] 2020 * obs obs 2021 */ 2022 for (i=0;i<numgrids;i++){ 2023 misfit_list[i]=0.5*pow(meanvel,(double)2)*( 2024 pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) + 2025 pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) ); 2026 } 2027 2028 /*Process units: */ 2029 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,MisfitEnum,this->parameters); 2030 2031 } 2032 else{ 2033 /*Not supported yet! : */ 2034 ISSMERROR("%s%g","unsupported type of fit: ",fit); 2035 } 1933 1934 /*We are using an absolute misfit: 1935 * 1936 * 1 [ 2 2 ] 1937 * J = --- | (u - u ) + (v - v ) | 1938 * 2 [ obs obs ] 1939 * 1940 */ 1941 for (i=0;i<numgrids;i++){ 1942 misfit_list[i]=0.5*(pow((vx_list[i]-obs_vx_list[i]),(double)2)+pow((vy_list[i]-obs_vy_list[i]),(double)2)); 1943 } 1944 /*Process units: */ 1945 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceAbsVelMisfitEnum,this->parameters); 2036 1946 2037 1947 /*Apply weights to misfits*/ … … 2057 1967 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3); 2058 1968 2059 /*compute Misfit*/1969 /*compute SurfaceAbsVelMisfit*/ 2060 1970 Jelem+=misfit*Jdet*gauss_weight; 2061 1971 } … … 2066 1976 xfree((void**)&gauss_weights); 2067 1977 1978 /*Return: */ 1979 return Jelem; 1980 } 1981 /*}}}*/ 1982 /*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ 1983 double Tria::SurfaceRelVelMisfit(bool process_units){ 1984 1985 int i; 1986 1987 /* output: */ 1988 double Jelem=0; 1989 1990 /* node data: */ 1991 const int numgrids=3; 1992 const int numdof=2*numgrids; 1993 const int NDOF2=2; 1994 double xyz_list[numgrids][3]; 1995 1996 /* grid data: */ 1997 double vx_list[numgrids]; 1998 double vy_list[numgrids]; 1999 double obs_vx_list[numgrids]; 2000 double obs_vy_list[numgrids]; 2001 double misfit_square_list[numgrids]; 2002 double misfit_list[numgrids]; 2003 double weights_list[numgrids]; 2004 2005 /* gaussian points: */ 2006 int num_gauss,ig; 2007 double* first_gauss_area_coord = NULL; 2008 double* second_gauss_area_coord = NULL; 2009 double* third_gauss_area_coord = NULL; 2010 double* gauss_weights = NULL; 2011 double gauss_weight; 2012 double gauss_l1l2l3[3]; 2013 double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; 2014 2015 /* parameters: */ 2016 double velocity_mag,obs_velocity_mag; 2017 double misfit; 2018 2019 /* Jacobian: */ 2020 double Jdet; 2021 2022 /*relative and logarithmic control method :*/ 2023 double scalex=1; 2024 double scaley=1; 2025 2026 /*inputs: */ 2027 bool onwater; 2028 2029 double meanvel, epsvel; 2030 2031 /*retrieve inputs :*/ 2032 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2033 2034 /*If on water, return 0: */ 2035 if(onwater)return 0; 2036 2037 /* Get node coordinates and dof list: */ 2038 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2039 2040 /* Recover input data: */ 2041 inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); 2042 inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); 2043 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); 2044 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); 2045 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 2046 2047 /*retrieve some parameters: */ 2048 this->parameters->FindParam(&meanvel,MeanVelEnum); 2049 this->parameters->FindParam(&epsvel,EpsVelEnum); 2050 2051 /* Compute SurfaceRelVelMisfit at the 3 nodes 2052 * Here we integrate linearized functions: 2053 * 2054 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2055 * 2056 * where J_i are the misfits at the 3 nodes of the triangle 2057 * Phi_i is the nodal function (P1) with respect to 2058 * the vertex i 2059 */ 2060 2061 /*We are using a relative misfit: 2062 * 2063 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 2064 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 2065 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 2066 * obs obs 2067 */ 2068 for (i=0;i<numgrids;i++){ 2069 scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2); 2070 scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2); 2071 if(obs_vx_list[i]==0)scalex=0; 2072 if(obs_vy_list[i]==0)scaley=0; 2073 misfit_list[i]=0.5*(scalex*pow((vx_list[i]-obs_vx_list[i]),2)+scaley*pow((vy_list[i]-obs_vy_list[i]),2)); 2074 } 2075 2076 /*Process units: */ 2077 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceRelVelMisfitEnum,this->parameters); 2078 2079 /*Apply weights to misfits*/ 2080 for (i=0;i<numgrids;i++){ 2081 misfit_list[i]=weights_list[i]*misfit_list[i]; 2082 } 2083 2084 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 2085 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 2086 2087 /* Start looping on the number of gaussian points: */ 2088 for (ig=0; ig<num_gauss; ig++){ 2089 /*Pick up the gaussian point: */ 2090 gauss_weight=*(gauss_weights+ig); 2091 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 2092 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 2093 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 2094 2095 /* Get Jacobian determinant: */ 2096 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); 2097 2098 /*Compute misfit at gaussian point: */ 2099 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3); 2100 2101 /*compute SurfaceRelVelMisfit*/ 2102 Jelem+=misfit*Jdet*gauss_weight; 2103 } 2104 2105 xfree((void**)&first_gauss_area_coord); 2106 xfree((void**)&second_gauss_area_coord); 2107 xfree((void**)&third_gauss_area_coord); 2108 xfree((void**)&gauss_weights); 2109 2110 /*Return: */ 2111 return Jelem; 2112 } 2113 /*}}}*/ 2114 /*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/ 2115 double Tria::SurfaceLogVelMisfit(bool process_units){ 2116 2117 int i; 2118 2119 /* output: */ 2120 double Jelem=0; 2121 2122 /* node data: */ 2123 const int numgrids=3; 2124 const int numdof=2*numgrids; 2125 const int NDOF2=2; 2126 double xyz_list[numgrids][3]; 2127 2128 /* grid data: */ 2129 double vx_list[numgrids]; 2130 double vy_list[numgrids]; 2131 double obs_vx_list[numgrids]; 2132 double obs_vy_list[numgrids]; 2133 double misfit_square_list[numgrids]; 2134 double misfit_list[numgrids]; 2135 double weights_list[numgrids]; 2136 2137 /* gaussian points: */ 2138 int num_gauss,ig; 2139 double* first_gauss_area_coord = NULL; 2140 double* second_gauss_area_coord = NULL; 2141 double* third_gauss_area_coord = NULL; 2142 double* gauss_weights = NULL; 2143 double gauss_weight; 2144 double gauss_l1l2l3[3]; 2145 double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; 2146 2147 /* parameters: */ 2148 double velocity_mag,obs_velocity_mag; 2149 double misfit; 2150 2151 /* Jacobian: */ 2152 double Jdet; 2153 2154 /*relative and logarithmic control method :*/ 2155 double scalex=1; 2156 double scaley=1; 2157 2158 /*inputs: */ 2159 bool onwater; 2160 2161 double meanvel, epsvel; 2162 2163 /*retrieve inputs :*/ 2164 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2165 2166 /*If on water, return 0: */ 2167 if(onwater)return 0; 2168 2169 /* Get node coordinates and dof list: */ 2170 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2171 2172 /* Recover input data: */ 2173 inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); 2174 inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); 2175 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); 2176 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); 2177 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 2178 2179 /*retrieve some parameters: */ 2180 this->parameters->FindParam(&meanvel,MeanVelEnum); 2181 this->parameters->FindParam(&epsvel,EpsVelEnum); 2182 2183 /* Compute SurfaceLogVelMisfit at the 3 nodes 2184 * Here we integrate linearized functions: 2185 * 2186 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2187 * 2188 * where J_i are the misfits at the 3 nodes of the triangle 2189 * Phi_i is the nodal function (P1) with respect to 2190 * the vertex i 2191 */ 2192 2193 /*We are using a logarithmic misfit: 2194 * 2195 * [ vel + eps ] 2 2196 * J = 4 \bar{v}^2 | log ( ----------- ) | 2197 * [ vel + eps ] 2198 * obs 2199 */ 2200 for (i=0;i<numgrids;i++){ 2201 velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil. 2202 obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil. 2203 misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); 2204 } 2205 2206 /*Process units: */ 2207 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceLogVelMisfitEnum,this->parameters); 2208 2209 /*Apply weights to misfits*/ 2210 for (i=0;i<numgrids;i++){ 2211 misfit_list[i]=weights_list[i]*misfit_list[i]; 2212 } 2213 2214 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 2215 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 2216 2217 /* Start looping on the number of gaussian points: */ 2218 for (ig=0; ig<num_gauss; ig++){ 2219 /*Pick up the gaussian point: */ 2220 gauss_weight=*(gauss_weights+ig); 2221 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 2222 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 2223 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 2224 2225 /* Get Jacobian determinant: */ 2226 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); 2227 2228 /*Compute misfit at gaussian point: */ 2229 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3); 2230 2231 /*compute SurfaceLogVelMisfit*/ 2232 Jelem+=misfit*Jdet*gauss_weight; 2233 } 2234 2235 xfree((void**)&first_gauss_area_coord); 2236 xfree((void**)&second_gauss_area_coord); 2237 xfree((void**)&third_gauss_area_coord); 2238 xfree((void**)&gauss_weights); 2239 2240 /*Return: */ 2241 return Jelem; 2242 } 2243 /*}}}*/ 2244 /*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/ 2245 double Tria::SurfaceLogVxVyMisfit(bool process_units){ 2246 2247 int i; 2248 2249 /* output: */ 2250 double Jelem=0; 2251 2252 /* node data: */ 2253 const int numgrids=3; 2254 const int numdof=2*numgrids; 2255 const int NDOF2=2; 2256 double xyz_list[numgrids][3]; 2257 2258 /* grid data: */ 2259 double vx_list[numgrids]; 2260 double vy_list[numgrids]; 2261 double obs_vx_list[numgrids]; 2262 double obs_vy_list[numgrids]; 2263 double misfit_square_list[numgrids]; 2264 double misfit_list[numgrids]; 2265 double weights_list[numgrids]; 2266 2267 /* gaussian points: */ 2268 int num_gauss,ig; 2269 double* first_gauss_area_coord = NULL; 2270 double* second_gauss_area_coord = NULL; 2271 double* third_gauss_area_coord = NULL; 2272 double* gauss_weights = NULL; 2273 double gauss_weight; 2274 double gauss_l1l2l3[3]; 2275 double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; 2276 2277 /* parameters: */ 2278 double velocity_mag,obs_velocity_mag; 2279 double misfit; 2280 2281 /* Jacobian: */ 2282 double Jdet; 2283 2284 /*relative and logarithmic control method :*/ 2285 double scalex=1; 2286 double scaley=1; 2287 double S=0; 2288 int fit=-1; 2289 2290 /*inputs: */ 2291 bool onwater; 2292 2293 double meanvel, epsvel; 2294 2295 /*retrieve inputs :*/ 2296 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2297 2298 /*If on water, return 0: */ 2299 if(onwater)return 0; 2300 2301 /* Get node coordinates and dof list: */ 2302 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2303 2304 /* Recover input data: */ 2305 inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); 2306 inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); 2307 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); 2308 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); 2309 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 2310 2311 /*retrieve some parameters: */ 2312 this->parameters->FindParam(&meanvel,MeanVelEnum); 2313 this->parameters->FindParam(&epsvel,EpsVelEnum); 2314 2315 /* Compute SurfaceLogVxVyMisfit at the 3 nodes 2316 * Here we integrate linearized functions: 2317 * 2318 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2319 * 2320 * where J_i are the misfits at the 3 nodes of the triangle 2321 * Phi_i is the nodal function (P1) with respect to 2322 * the vertex i 2323 */ 2324 2325 /*We are using an logarithmic 2 misfit: 2326 * 2327 * 1 [ |u| + eps 2 |v| + eps 2 ] 2328 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 2329 * 2 [ |u |+ eps |v |+ eps ] 2330 * obs obs 2331 */ 2332 for (i=0;i<numgrids;i++){ 2333 misfit_list[i]=0.5*pow(meanvel,(double)2)*( 2334 pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) + 2335 pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) ); 2336 } 2337 2338 /*Process units: */ 2339 if(process_units)NodalValuesUnitConversion(&misfit_list[0],numgrids,SurfaceLogVxVyMisfitEnum,this->parameters); 2340 2341 /*Apply weights to misfits*/ 2342 for (i=0;i<numgrids;i++){ 2343 misfit_list[i]=weights_list[i]*misfit_list[i]; 2344 } 2345 2346 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 2347 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 2348 2349 /* Start looping on the number of gaussian points: */ 2350 for (ig=0; ig<num_gauss; ig++){ 2351 /*Pick up the gaussian point: */ 2352 gauss_weight=*(gauss_weights+ig); 2353 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 2354 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 2355 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 2356 2357 /* Get Jacobian determinant: */ 2358 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); 2359 2360 /*Compute misfit at gaussian point: */ 2361 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3); 2362 2363 /*compute SurfaceLogVxVyMisfit*/ 2364 Jelem+=misfit*Jdet*gauss_weight; 2365 } 2366 2367 xfree((void**)&first_gauss_area_coord); 2368 xfree((void**)&second_gauss_area_coord); 2369 xfree((void**)&third_gauss_area_coord); 2370 xfree((void**)&gauss_weights); 2371 2372 /*Return: */ 2373 return Jelem; 2374 } 2375 /*}}}*/ 2376 /*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/ 2377 double Tria::SurfaceAverageVelMisfit(bool process_units){ 2378 2379 int i; 2380 2381 /* output: */ 2382 double Jelem=0; 2383 2384 /* node data: */ 2385 const int numgrids=3; 2386 const int numdof=2*numgrids; 2387 const int NDOF2=2; 2388 double xyz_list[numgrids][3]; 2389 2390 /* grid data: */ 2391 double vx_list[numgrids]; 2392 double vy_list[numgrids]; 2393 double obs_vx_list[numgrids]; 2394 double obs_vy_list[numgrids]; 2395 double misfit_square_list[numgrids]; 2396 double misfit_list[numgrids]; 2397 double weights_list[numgrids]; 2398 2399 /* gaussian points: */ 2400 int num_gauss,ig; 2401 double* first_gauss_area_coord = NULL; 2402 double* second_gauss_area_coord = NULL; 2403 double* third_gauss_area_coord = NULL; 2404 double* gauss_weights = NULL; 2405 double gauss_weight; 2406 double gauss_l1l2l3[3]; 2407 double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; 2408 2409 /* parameters: */ 2410 double velocity_mag,obs_velocity_mag; 2411 double misfit; 2412 2413 /* Jacobian: */ 2414 double Jdet; 2415 2416 /*relative and logarithmic control method :*/ 2417 double scalex=1; 2418 double scaley=1; 2419 double S=0; 2420 int fit=-1; 2421 2422 /*inputs: */ 2423 bool onwater; 2424 2425 double meanvel, epsvel; 2426 2427 /*retrieve inputs :*/ 2428 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2429 2430 /*If on water, return 0: */ 2431 if(onwater)return 0; 2432 2433 /* Get node coordinates and dof list: */ 2434 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); 2435 2436 /* Recover input data: */ 2437 inputs->GetParameterValue(&S,SurfaceAreaEnum); 2438 inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); 2439 inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); 2440 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); 2441 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); 2442 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 2443 2444 /*retrieve some parameters: */ 2445 this->parameters->FindParam(&meanvel,MeanVelEnum); 2446 this->parameters->FindParam(&epsvel,EpsVelEnum); 2447 2448 /* Compute SurfaceAverageVelMisfit at the 3 nodes 2449 * Here we integrate linearized functions: 2450 * 2451 * J(E) = int_E sum_{i=1}^3 J_i Phi_i 2452 * 2453 * where J_i are the misfits at the 3 nodes of the triangle 2454 * Phi_i is the nodal function (P1) with respect to 2455 * the vertex i 2456 */ 2457 2458 /*We are using a spacially average absolute misfit: 2459 * 2460 * 1 2 2 2461 * J = --- sqrt( (u - u ) + (v - v ) ) 2462 * S obs obs 2463 */ 2464 for (i=0;i<numgrids;i++) misfit_square_list[i]=pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2); 2465 2466 /*Process units: */ 2467 if(process_units)NodalValuesUnitConversion(&misfit_square_list[0],numgrids,SurfaceAverageVelMisfitEnum,this->parameters); 2468 2469 /*Take the square root, and scale by surface: */ 2470 for (i=0;i<numgrids;i++)misfit_list[i]=pow(misfit_square_list[i],2)/S; 2471 2472 /*Apply weights to misfits*/ 2473 for (i=0;i<numgrids;i++){ 2474 misfit_list[i]=weights_list[i]*misfit_list[i]; 2475 } 2476 2477 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 2478 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); 2479 2480 /* Start looping on the number of gaussian points: */ 2481 for (ig=0; ig<num_gauss; ig++){ 2482 /*Pick up the gaussian point: */ 2483 gauss_weight=*(gauss_weights+ig); 2484 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 2485 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig); 2486 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 2487 2488 /* Get Jacobian determinant: */ 2489 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); 2490 2491 /*Compute misfit at gaussian point: */ 2492 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3); 2493 2494 /*compute SurfaceAverageVelMisfit*/ 2495 Jelem+=misfit*Jdet*gauss_weight; 2496 } 2497 2498 xfree((void**)&first_gauss_area_coord); 2499 xfree((void**)&second_gauss_area_coord); 2500 xfree((void**)&third_gauss_area_coord); 2501 xfree((void**)&gauss_weights); 2502 2068 2503 /*Return: */ 2069 2504 return Jelem; … … 2156 2591 /*inputs: */ 2157 2592 bool onwater; 2158 int fit;2159 2593 2160 2594 /*retrieve inputs :*/ 2161 inputs->GetParameterValue(&fit,FitEnum);2162 2595 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2163 2164 /*If fit!=3, do not compute surface: */2165 if(fit!=3)return 0;2166 2596 2167 2597 /*If on water, return 0: */ … … 4471 4901 double scale=0; 4472 4902 double S=0; 4473 int fit=-1;4903 int response; 4474 4904 4475 4905 /*retrieve some parameters: */ … … 4490 4920 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 4491 4921 4492 inputs->GetParameterValue(& fit,FitEnum);4493 if( fit==3){4922 inputs->GetParameterValue(&response,CmResponseEnum); 4923 if(response==SurfaceAverageVelMisfitEnum){ 4494 4924 inputs->GetParameterValue(&S,SurfaceAreaEnum); 4495 4925 } … … 4508 4938 * the vertex i 4509 4939 */ 4510 if( fit==0){4940 if(response==SurfaceAbsVelMisfitEnum){ 4511 4941 /*We are using an absolute misfit: 4512 4942 * … … 4524 4954 } 4525 4955 } 4526 else if( fit==1){4956 else if(response==SurfaceRelVelMisfitEnum){ 4527 4957 /*We are using a relative misfit: 4528 4958 * … … 4546 4976 } 4547 4977 } 4548 else if( fit==2){4978 else if(response==SurfaceLogVelMisfitEnum){ 4549 4979 /*We are using a logarithmic misfit: 4550 4980 * … … 4567 4997 } 4568 4998 } 4569 else if( fit==3){4999 else if(response==SurfaceAverageVelMisfitEnum){ 4570 5000 /*We are using an spacially average absolute misfit: 4571 5001 * … … 4584 5014 } 4585 5015 } 4586 else if( fit==4){5016 else if(response==SurfaceLogVxVyMisfitEnum){ 4587 5017 /*We are using an logarithmic 2 misfit: 4588 5018 * … … 4604 5034 else{ 4605 5035 /*Not supported yet! : */ 4606 ISSMERROR(" %s%g","unsupported type of fit: ",fit);5036 ISSMERROR("response %s not supported yet",EnumToString(response)); 4607 5037 } 4608 5038 … … 4710 5140 double S=0; 4711 5141 int fit=-1; 5142 int response; 4712 5143 4713 5144 /*retrieve some parameters: */ … … 4728 5159 inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); 4729 5160 4730 inputs->GetParameterValue(& fit,FitEnum);4731 if( fit==3){5161 inputs->GetParameterValue(&response,CmResponseEnum); 5162 if(response==SurfaceAverageVelMisfitEnum){ 4732 5163 inputs->GetParameterValue(&S,SurfaceAreaEnum); 4733 5164 } … … 4746 5177 * the vertex i 4747 5178 */ 4748 if( fit==0){5179 if(response==SurfaceAbsVelMisfitEnum){ 4749 5180 /*We are using an absolute misfit: 4750 5181 * … … 4762 5193 } 4763 5194 } 4764 else if( fit==1){5195 else if(response==SurfaceRelVelMisfitEnum){ 4765 5196 /*We are using a relative misfit: 4766 5197 * … … 4784 5215 } 4785 5216 } 4786 else if( fit==2){5217 else if(response==SurfaceLogVelMisfitEnum){ 4787 5218 /*We are using a logarithmic misfit: 4788 5219 * … … 4805 5236 } 4806 5237 } 4807 else if( fit==3){5238 else if(response==SurfaceAverageVelMisfitEnum){ 4808 5239 /*We are using an spacially average absolute misfit: 4809 5240 * … … 4822 5253 } 4823 5254 } 4824 else if( fit==4){5255 else if(response==SurfaceLogVxVyMisfitEnum){ 4825 5256 /*We are using an logarithmic 2 misfit: 4826 5257 * … … 4842 5273 else{ 4843 5274 /*Not supported yet! : */ 4844 ISSMERROR(" %s%g","unsupported type of fit: ",fit);5275 ISSMERROR("response %s not supported yet",EnumToString(response)); 4845 5276 } 4846 5277 … … 6298 6729 name==VxObsEnum || 6299 6730 name==VyObsEnum || 6300 name== FitEnum ||6731 name==CmResponseEnum || 6301 6732 name==DragCoefficientEnum || 6302 6733 name==GradientEnum || -
issm/trunk/src/c/objects/Elements/Tria.h
r5258 r5281 68 68 void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); 69 69 void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); 70 double CostFunction(bool process_units);70 double RegularizeInversion(void); 71 71 void CreateKMatrix(Mat Kgg); 72 72 void CreatePVector(Vec pg); … … 103 103 void MinVy(double* pminvy, bool process_units); 104 104 void MinVz(double* pminvz, bool process_units); 105 double Misfit(bool process_units); 105 double SurfaceAbsVelMisfit(bool process_units); 106 double SurfaceRelVelMisfit(bool process_units); 107 double SurfaceLogVelMisfit(bool process_units); 108 double SurfaceLogVxVyMisfit(bool process_units); 109 double SurfaceAverageVelMisfit(bool process_units); 106 110 void PatchFill(int* pcount, Patch* patch); 107 111 void PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes); -
issm/trunk/src/c/objects/IoModel.cpp
r5251 r5281 86 86 xfree((void**)&this->rheology_B); 87 87 xfree((void**)&this->rheology_n); 88 xfree((void**)&this-> fit);88 xfree((void**)&this->cm_responses); 89 89 xfree((void**)&this->weights); 90 90 xfree((void**)&this->cm_jump); … … 302 302 303 303 /*!solution parameters: */ 304 this-> fit=NULL;304 this->cm_responses=NULL; 305 305 this->weights=NULL; 306 306 this->cm_jump=NULL; -
issm/trunk/src/c/objects/IoModel.h
r5251 r5281 129 129 130 130 /*solution parameters: */ 131 double* fit;131 double* cm_responses; 132 132 double* weights; 133 133 double* cm_jump; -
issm/trunk/src/c/solutions/control_core.cpp
r5247 r5281 26 26 int solution_type; 27 27 28 double* fit=NULL;28 double* responses=NULL; 29 29 double* optscal=NULL; 30 30 double* maxiter=NULL; … … 46 46 femmodel->parameters->FindParam(&nsteps,NStepsEnum); 47 47 femmodel->parameters->FindParam(&control_type,ControlTypeEnum); 48 femmodel->parameters->FindParam(& fit,NULL,FitEnum);48 femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); 49 49 femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); 50 50 femmodel->parameters->FindParam(&maxiter,NULL,MaxIterEnum); … … 65 65 } 66 66 67 /*Initialize misfit: */67 /*Initialize responses: */ 68 68 J=(double*)xmalloc(nsteps*sizeof(double)); 69 69 … … 75 75 for(n=0;n<nsteps;n++){ 76 76 77 /*Display info*/ 77 78 _printf_("\n%s%i%s%i\n"," control method step ",n+1,"/",nsteps); 78 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int) fit[n],FitEnum);79 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n],CmResponseEnum); 79 80 80 81 /*In case we are running a steady state control method, compute new temperature field using new parameter distribution: */ … … 97 98 InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar*optscal[n],true); 98 99 99 if(controlconvergence(J, fit,eps_cm,n)) break;100 if(controlconvergence(J,responses,eps_cm,n)) break; 100 101 101 102 /*Temporary saving every 5 control steps: */ … … 117 118 cleanup_and_return: 118 119 /*Free ressources: */ 119 xfree((void**)& fit);120 xfree((void**)&responses); 120 121 xfree((void**)&optscal); 121 122 xfree((void**)&maxiter); -
issm/trunk/src/c/solutions/controlconvergence.cpp
r4778 r5281 18 18 19 19 20 bool controlconvergence(double* J, double* fit, double eps_cm, int n){20 bool controlconvergence(double* J,double* responses, double eps_cm, int n){ 21 21 22 22 int i; … … 29 29 //go through the previous misfits(starting from n-2) 30 30 while(i>=0){ 31 if ( fit[i]==fit[n]){31 if (responses[i]==responses[n]){ 32 32 //convergence test only if we have the same misfits 33 33 if ((J[i]-J[n])/J[n] <= eps_cm){ -
issm/trunk/src/c/solutions/objectivefunctionC.cpp
r5251 r5281 28 28 29 29 /*parameters: */ 30 FemModel *femmodel = NULL; 31 int n; 32 double *optscal = NULL; 33 double *fit = NULL; 34 int control_type; 35 int solution_type; 36 int analysis_type; 37 bool isstokes = false; 38 bool conserve_loads = true; 39 bool process_units = false; 40 30 FemModel *femmodel = NULL; 31 int n; 32 double *optscal = NULL; 33 double *responses =NULL; 34 int control_type; 35 int solution_type; 36 int analysis_type; 37 bool isstokes = false; 38 bool conserve_loads = true; 39 41 40 /*Recover finite element model: */ 42 41 femmodel=optargs->femmodel; … … 46 45 47 46 femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); 48 femmodel->parameters->FindParam(& fit,NULL,FitEnum);47 femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); 49 48 femmodel->parameters->FindParam(&isstokes,IsStokesEnum); 50 49 femmodel->parameters->FindParam(&control_type,ControlTypeEnum); … … 82 81 83 82 /*Compute misfit for this velocity field.*/ 84 InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,(int)fit[n],FitEnum); 85 CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,process_units); 83 CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,(int)responses[n]); 86 84 87 85 /*Free ressources:*/ 88 xfree((void**)& fit);86 xfree((void**)&responses); 89 87 xfree((void**)&optscal); 90 88 return J; -
issm/trunk/src/c/solutions/solutions.h
r5247 r5281 32 32 //convergence: 33 33 void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters); 34 bool controlconvergence(double* J, double* fit, double eps_cm, int n);34 bool controlconvergence(double* J, double* response, double eps_cm, int n); 35 35 bool steadystateconvergence(FemModel* femmodel); 36 36
Note:
See TracChangeset
for help on using the changeset viewer.
